0

I am solving an optimization problem using MATLAB with YALMIP which is similar to a unit commitment problem.

This is a Mixed Integer Linear Programming problem.

I am considering two decision variables --> P and T_room

I formulated the problem in MATLAB with YALMIP,

When I use CPLEX and MOSEK solvers, I am getting different results.

Can anyone explain to me why that is happening?

CPLEX solver results

MOSEK solver results

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([6,8],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);

T_room_out = value(T_room);
  • 1
    Are the results just different in detail, but 'equivalent' quality, i.e. do they have the same objective value? Very often there are many equivalent solutions which corespond to symmetries and or different permutations of the speciific choices. – TimChippingtonDerrick Dec 06 '19 at 10:57
  • 1
    In case the results are qualitatively different (different objective value) then the first thing to try is to export both models to an LP file and see if they are the same. – Daniel Junglas Dec 09 '19 at 08:17
  • @TimChippingtonDerrick Thank you for the response. The thing is usually MOSEK takes a long time to solve the problem compared to CPLEX. And regarding the quality of the output, both of them gives different results. – gayan_lanke Dec 09 '19 at 23:45
  • @DanielJunglas Thank you. Can you just provide more details regarding this. – gayan_lanke Dec 09 '19 at 23:46
  • Sorry, I don't know YALMIP. You may want to ask a YALMIP expert about this. Note that you still not answer what are the differences in the solutions. Do they really give different objective values? Or is it just different solution vectors that have the same objective value. – Daniel Junglas Dec 10 '19 at 12:45

1 Answers1

1

Your code does not make sense to begin with, as your objective is a vector

>> Objective
Linear matrix variable 3x1 (full, real, 42 variables)

With that, YALMIP will solve three optimization problems, and then values you obtain when applying value command will be the values from the last solution.

Having said that, concentrating on the last objective you have defined (although I think you intended to minimize the sum or something), I get the same objective with both solvers (not the same solution, the solution simply isn't unique)

>> solution  = optimize (Constraints, Objective(3), sdpsettings('solver','cplex'))
>> value(Objective(3))
ans =
    1.8418
>> value(P)
ans =
         0    5.2500    3.5000    5.2500    3.5000    5.2500    5.2500    5.2500
         0    3.0000    3.0000    6.0000    6.0000    6.0000    3.0000    4.5000
    5.0000    5.0000    5.0000    5.0000    5.0000    5.0000         0         0
>> solution  = optimize (Constraints, Objective(3), sdpsettings('solver','mosek'))
>> value(Objective(3))
ans =
    1.8418
ans =
    3.5000    3.5000    3.5000    3.5000    3.5000    3.5000    5.2500    5.2500
    3.0000    3.0000    3.0000    3.0000    6.0000         0    3.0000    6.0000
    3.7500    5.0000    5.0000    5.0000    5.0000    5.0000         0         0
Johan Löfberg
  • 378
  • 1
  • 13