Regarding infeasibility of a problem solved with cplex

[Gayan Lankeshwara]

Hi !

I encountered the following issue for an optimization problem (similar to unit commitment) solved with cplex.

And found that it is something to deal with my temperature model constraints.

Temeprature model is attached herewith.

There I know all the other parameters except T_room,i,t  .

At the same time, the intial vector of T_room (for the 1st time step) is also avaialable with me.

Can someone please tell me how to form the constraints for T_room considering the above ?

Thanks.

[Johan Löfberg]

How is it any different from all other MPC models etc in the Wiki?

Troom(i,t) == Tout(i,t)-Q(i,t)*R-...

[Gayan Lankeshwara]

Hi Johan,

I am kind of struggling to solve the attached optimisation problem. It says "Infeasible problem (CPLEX-IBM)"

Once I remove that constraint, it works but does not make sense.

I will attach the code as well as the problem formulation as well.

Hope I can get a help from you.

Thanks a lot.

[Johan Löfberg]

It's simply not feasible. Easily seen by e.g. adding slacks on your T_room constraints and the minimal possible value is around 5.4

[Gayan Lankeshwara]

Hi Johan !

Can you just explain this a little bit.

I know that there is something wrong with the T_room formulation. But I did not get it right ?

Thanks.

[Johan Löfberg]

There's simply no solution to the problem. Heck, already this is infeasible

Constraints = [];

for k = 1: Horizon    

    Constraints = [Constraints,  P(:,k) <= Pmax];

end

for k = 2: Horizon    

    Constraints = [Constraints, T_lower<=T_room(:,k)<=T_upper];

end

for k = 2: Horizon    

  Constraints = [Constraints, T_room(:,k)+(c/a).*(1 - exp(-(deltaT./(R.*Ca)))).*P(:,k)-T_room(:,k-1).*exp(-deltaT./(R.*Ca)) == T_outdoor - ((a*c - b*d)/a).*R.*(1 - exp(-(deltaT./(R.*Ca))))-T_outdoor.* exp(-deltaT./(R.*Ca))];    

end

optimize(Constraints)

You don't have enough cooling power to keep all rooms within specified limits, regardless of initial condition in the rooms!

Adding a slack on the available power shows that the smallest possible slack is 4.57, i.e. you need roughly twice as much power

[Gayan Lankeshwara]

Hi Johan,

Thank you for the response.

I understood what you said, but still I am struggling with the same problem.

I updated the constraint on room temperature and now the problem is feasible.

But I have the following issues.

1.  I developed this problem with the help of Unit commitment problem in YALMIP. There I see the objective function is minimizing the total cost in the horizon with the following .

Objective = Objective + P(:,k)'*Q*P(:,k) + C*P(:,k);

Is this the same as minimizing the cost in each time interval (delta T ) or else, is there a difference between minimizing the cost in the total horizon vs, minimizing the cost in each time interval.

If there is a difference, how can I change the optimization problem to minimize the cost in each time interval 

2. I developed the objective function with a linear cost and a linear model for the inconvenience as in the code attached (with here). When I run with both CPLEX and MOSEK, they give different results which I could not understand why?

   At the same time, both the solutions seem to violate the following condition.

for k = 1:Horizon

  Constraints = [Constraints, sum(P(:,k)) >= Pdemand(k)];
end

I tried with a Demand slack variable with penalty levels, but once i do it, the problem turns to be infeasible.

3. Then I tried solving the unit commitment example with the same quadratic function but with quantised power levels for all the 3 generators.

 But it seems to take a long time to solve in my desktop. Why is this ?

Nunits = 3;
Horizon = 48;


Pmax = [100;50;25];
Pmin = [20;40;1];


Unit1Levels = [0  50 75 100];
Unit2Levels = [0 10 20 30 40 50];
Unit3Levels = [0 1 6 10 12 20];


Q = diag([.04 .01 .02]);
C = [10 20 20]


Pforecast = 100 + 50*sin((1:Horizon)*2*pi/24)


onoff = binvar(Nunits,Horizon,'full');
P     = sdpvar(Nunits,Horizon,'full');

Constraints = [];


for k = 1:Horizon

  Constraints = [Constraints, onoff(:,k).*Pmin <= P(:,k) <= onoff(:,k).*Pmax];
end

for k = 1:Horizon

  Constraints = [Constraints, sum(P(:,k)) >= Pforecast(k)];
end


%% different quantised levels for the appliances

for k = 1:Horizon

  Constraints = [Constraints, ismember(P(1,k),Unit1Levels)];
end

for k = 1:Horizon

      Constraints = [Constraints, ismember(P(2,k),Unit2Levels)];
end

for k = 1:Horizon

      Constraints = [Constraints, ismember(P(3,k),Unit3Levels)];
end


Objective = 0;

for k = 1:Horizon

  Objective = Objective + P(:,k)'*Q*P(:,k) + C*P(:,k);
end


ops = sdpsettings('verbose',2,'debug',1, 'solver', 'mosek');
optimize(Constraints,Objective,ops)


value(P)
value(onoff)

I would really appreciate if you could help me in the above.

Thanks in advance,

BR,

Gayan

[Gayan Lankeshwara]

I missed to attach the code for Q2.

Here is it.

%%% 03.12.2019




%{


Here we consider the approach where the two objectives are taken into
consideration.


For inverter type air condtioners with different power level reductions


One big asumption is that all the air condtioners are of the same type (same manufacturer and same type)


We consider the same temperature for all the iterations which needs to be
changed in the next version




Note: The issue here is the problem is infeasible.


Need to think how to resolve it.


%% Assuming a very small distribution of air conditioners N = 2




%}
clc;
clear all;


rng(0,'twister');




%% time step of the demand response --> 15 mins


deltaT = 15/60 ;


%% number of air conditioners
Nunits = 3;


%% demand required for each interval


Pdemand = [2 3 4 2 5 5 8 9]';


% Pdemand = randi([2,10],8,1);


%% time interval considered for the study --> 15 mins intervals for 2
%%% hours


Horizon = length(Pdemand) ;




%% cost for each appliance ( need to think little about this approach )


Q = diag([rand(1,Nunits)]);


%% defining the possible reduction levels for the appliances


% ratedPower  = [1 6 4]';


ratedPower  = randi([5,7],Nunits,1);


%% defining minimum and maximum power limits for the appliances


Pmax = randi([8,9],Nunits,1);
Pmin = zeros(Nunits,1);


% possible power levels for  appliances (discrete levels)


LevelsOutput = [0 0.5 0.75 1.0] ;


% matrix containing all the possible levels of output power for all the
%% appliances


PowerLevelMatrix = ratedPower * LevelsOutput ;


%% parameters based on the temperature


%%% outside temperature (we assume that the outdoor temperature remains the same throughout the optimisation process)


T_outdoor = randi([34, 38], Nunits,1);


%%% preferred lower level of temperature for the appliances
T_lower  =  randi([20,21], Nunits,1) ;


%%% preferred upper level of temperature for the appliances
T_upper  =  randi([28, 29], Nunits,1) ;


%%% initial temperature for the appliances
T_initial = randi([22, 28], Nunits,1);




%% parameters of the house


%% equivalent thermal resistance
R = (2.5-1.5)*rand(Nunits,1) + 1.5 ;


%%% equivalent thermal capacitance
Ca = (2.5-1.5)*rand(Nunits,1) + 1.5 ;


%% for the air condtioner


c = repmat(0.03,Nunits,1);
d = repmat(-0.4,Nunits,1);
a = repmat(0.06,Nunits,1);
b = repmat(-0.3,Nunits,1);


%% defining the variables of the demand response problem


%%% power levels for the appliances


P = sdpvar(Nunits, Horizon, 'full');


%% slack variable for exceeding the demand


% Demand_Penalty  = 50;
% DemandSlack = sdpvar(1,Horizon);


%% decision variable for the room temperature


T_room =  sdpvar(Nunits,Horizon,'full') ;




%% Constraints for the optimisation problem

Constraints = [];


for k = 1: Horizon
   

%     Constraints = [Constraints, sum(P(:,k)) + DemandSlack(k) >= Pdemand(k)];
%     Constraints = [Constraints,DemandSlack(k) >= 0];

      Constraints = [Constraints, sum(P(:,k))>= Pdemand(k)];
end

%%% constraint on min amd max power levels

for k = 1: Horizon

   
    Constraints = [Constraints, Pmin <= P(:,k) <= Pmax];
end


%%% adding the constraint on discrete power levels for each appliance
for nApp =  1 : Nunits

   
    for k = 1 : Horizon

       
    Constraints = [Constraints, ismember(P(nApp,k),PowerLevelMatrix(nApp,:))];
       
       
    end
       
end




%% the temperature should be within the desired limits --> cannot go beyond those limits

for k = 2: Horizon
   
    Constraints = [Constraints, T_lower <= T_room(:,k) <= T_upper ];
   
end

%% adding the initial constraint of temperature


Constraints = [Constraints, T_room(:,1) == T_initial];




%%% adding the temperature model of the constraint for the inverter type
%%% air condtioner

for k = 2: Horizon

   
  Constraints = [Constraints, T_room(:,k) == T_outdoor - ((a./c).*P(:,k) + ((b.*c - a.*d)./c)).*R.*(1- exp((-deltaT)./(R.*Ca))) - T_outdoor.* exp((-deltaT)./(R.*Ca)) + T_room(:,k-1).* exp((-deltaT)./(R.*Ca)) ];
 
 
end




%% adding a sample constraints to the temperature update


% for k =2: Horizon
%    
%     Constraints = [Constraints, T_room(:,k) == T_room(:,k-1) ];
% end


%%% defining the objective function




% weighting factors


alpha = 1 ;
beta =  1 ;


Objective = 0 ;


for k = 2 : Horizon
   
      Objective  =  Objective + alpha * Q * P(:,k) + beta * ((2* T_room(:,k) - T_lower - T_upper)./(T_upper - T_lower));
   
%       Objective  =  Objective + P(:,k)' * Q * P(:,k);
     
%       Objective =   Objective +  Demand_Penalty * DemandSlack(k) ;
     
end




%%% solving the optimisation problem
%
ops = sdpsettings('verbose',2, 'debug', 1, 'solver', 'mosek');


% ops = sdpsettings('verbose',1, 'debug', 1);


solution  = optimize (Constraints, Objective, ops) ;




% Pout = value(P);


disp('The results with MOSEK solver')


disp('The output of P');


disp(value(P))


disp('The output of T_room');


T_room_out = value(T_room);


disp(value(T_room))


% stairs(Pout')
% ylim([0 10])

Thanks,

[Johan Löfberg]

1. What would it even mean to minimize the cost in every interval?? An objective is a scalar, and somehow all intervals have to boil down to a single measure. Whether you sum the squares, the absolute values, or take the largest value along the horizon is just a matter of norm

2. Define different....Have they solved the problem without any error or warning? Is the optimal objective the same 

3 https://yalmip.github.io/slowmilp