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?
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);