Problem Based approach for Mixed Integer Linear Programming problem fails to produce any feasible solution.

I am trying to implement a IEEE research paper for energy management in a Smart Home using MILP's problem based approach. I am doing the coding on the track of a similar example provided on mathworks site.

The link for that example is:- https://in.mathworks.com/videos/mixed-integer-linear-programming-in-matlab-91541.html?s_tid=srchtitle

The problem formulation for my problem is like this:-

The matlab coding that I have done is as follows:-

% Program for optimization and energy management in a Smart Home.
% Smart Home consists of a PV module, a wind turbine, a home battery...
% storage system, a Plug-in Electric Vehicle and an energy management system.
clc
% Load all the datas
LoadDemand=[0.1 0.5 0.4 3 5 3 3.5 0.1 0.1 1 3 2 2.3 2.3 2 2.5 5 6 7 6.5...
    4 2.5 2.3 0.1]; % Home Load Data for 24 hrs in kW
SolarIrradiation=[0 0 0 0 0.1 0.3 0.5 0.65 0.8 0.9 1 0.95 0.9 0.8 0.65...
    0.4 0.3 0.08 0 0 0 0 0 0]; % for 24 hrs in kW/m^2
WindVelocity=repmat(5,1,24); % for 24 hrs in m/sec
Cgrid=[0.033 0.027 0.020 0.017 0.017 0.029 0.033 0.054 0.215...
    0.572 0.572 0.572 0.215 0.572 0.286 0.279 0.086 0.059 0.050 0.061...
    0.181 0.077 0.043 0.037]; % all prices and costs are in Euros
OPcostPV=0.01; OPcostW=0.01; OPcostEVch=0.01; OPcostEVdisch=0.01;
OPcostBch=0.01; OPcostBdisch=0.01; SellingCost=0.29; % in Euros
% PreDefined variables
dt=1; % Time step
NOMb=10; % Nominal Battery Capacity in kWh
NOMbInit=6; % Initial Nominal Battery Capacity in kWh
SOCminB=0.2; % Battery State of Charge (minimum)
NOMev=24; % Nominal EV Battery Capacity in kWh
NOMevInit=16; % Initial Nominal EV Battery Capacity in kWh
SOCminev=0.2 % EV Battery State of Charge (minimum)
ec=1.1; % Battery Charging Coefficient
ed=0.9; % Battery Discharging Coefficient
PVeffi=0.186; % PV efficiency
A=25; % PV system area 
Vci=4; % Cut-in Speed for wind turbine in m/sec
Vco=24; % Cut-out Speed for wind turbine in m/sec
Vr=14; % rated speed for wind turbine in m/sec
PbChmax=1; % Battery max allowed charging power in kW
PbDischmax=1; % Battery max allowed discharging power in kW
PevChmax=3.3; % EV Battery max allowed charging power in kW
PevDischmax=3.3; % EV Battery max allowed charging power in kW
DevDrive=2; % Driving electricity demand of EV in kWh
% Initiate time variables for 24 hrs.
nHours=numel(LoadDemand);
Time=(1:nHours)';
idxHr2Toend=2:nHours;
DriveTime=Time(9:18); % Time for which ev is not at home
% Initiate the problem
powerprob=optimproblem;
% Initiate the variables:-
% Initaiate the continous variables
Pgrid=optimvar('Pgrid',nHours,1,'LowerBound',0,'UpperBound',5); % for grid power
Ppv=optimvar('Ppv',nHours,1,'LowerBound',0,'UpperBound',3.5); % for PV
Pwind=optimvar('Pwind',nHours,1,'LowerBound',0,'UpperBound',2.4); 
PevCh=optimvar('PevCh',nHours,1,'LowerBound',0,'UpperBound',3.3); % for ev charging
PevDisch=optimvar('PevDisch',nHours,1,'LowerBound',0,'UpperBound',3.3); % for ev discharge
PbCh=optimvar('PbCh',nHours,1,'LowerBound',0,'UpperBound',1); % for battery charge
PbDisch=optimvar('PbDisch',nHours,1,'LowerBound',0,'UpperBound',1); % for battery discharge
Psell=optimvar('Psell',nHours,1,'LowerBound',0,'UpperBound',inf); % for selling power to grid
SOCb=optimvar('SOCb',nHours,1,'LowerBound',0.2,'UpperBound',1); % Battery state of charge
SOCev=optimvar('SOCev',nHours,1,'LowerBound',0.2,'UpperBound',1); % EV Battery state of charge
% Initaiting the discrete variables
isONevCh=optimvar('isONevCh',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); % ON/OFF status of ev charge
isONevDisch=optimvar('isONevDisch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); % ON/OFF status of ev discharge
isONbCh=optimvar('isONbCh',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); % ON/OFF status of battery charge
isONbDisch=optimvar('isONbDisch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1); % ON/OFF status of battery charge
% Indicator if the home battery is starting to charge/discharge during the hour
startupBdisch=optimvar('startupBdisch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);
startupBch=optimvar('startupBch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);
% Indicator if the PEV battery is starting to charge/discharge during the hour
startupPevdisCh=optimvar('startupPEVdisch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);
startupPevCh=optimvar('startupPEVch',nHours,1,'Type','integer','LowerBound',0,'UpperBound',1);
% Defining the objective function cost variables
 gridCost = sum(((Pgrid*dt).*Cgrid'),1);

   PVGenCost = sum(((Ppv*dt)*OPcostPV),1);

   WindgenCost = sum(((Pwind*dt)*OPcostW),1);

   EVdischCost = sum(((PevDisch*dt)*OPcostEVdisch),1);

   BatterydischCost = sum(((PbDisch*dt)*OPcostBdisch),1);

   SellCost = sum(((Psell*dt)*SellingCost),1);
   % Set the objective function
  powerprob.Objective= gridCost + PVGenCost + WindgenCost + EVdischCost + BatterydischCost - SellCost;
  % All Equality Constraints :-
  % Power Balance Equality Constraint
  powerprob.Constraints.PowerBalance = Pgrid + Ppv + Pwind + PevDisch + PbDisch - LoadDemand' - PevCh - PbCh - Psell==0; 
  % Power Balance equation
  % Battery Stored Electricity Equality Constraint
  powerprob.Constraints.BatStoredElectricity=(NOMb*SOCb(idxHr2Toend,:))-(NOMb*SOCb(idxHr2Toend-1,:))-(((PbCh(idxHr2Toend,:)*dt)/ec)-(ed*PbDisch(idxHr2Toend,:)*dt))==0;
  powerprob.Constraints.InitBatStoredElectricity=(NOMb*SOCb(1))-NOMbInit-(((PbCh(idxHr2Toend,:)*dt)/ec)-(ed*PbDisch(idxHr2Toend,:)*dt))==0;
  % EV Battery Stored Electricity Equality Constraint
  powerprob.Constraints.EVBatStoredElectricity=(NOMev*SOCev(idxHr2Toend,:))-(NOMev*SOCev(idxHr2Toend-1,:))-(((PevCh(idxHr2Toend,:)*dt)/ec)-(ed*PevDisch(idxHr2Toend,:)*dt))==0;
  powerprob.Constraints.EVInitBatStoredElectricity=(NOMev*SOCev(1))-NOMevInit-(((PevCh(idxHr2Toend,:)*dt)/ec)-(ed*PevDisch(idxHr2Toend,:)*dt))==0;
  % Wind Generation Constraint
  if WindVelocity(1:nHours)<Vci 
   Pwind=0
  elseif WindVelocity(1:nHours)>Vco
      Pwind=0
  elseif Vr<WindVelocity(1:nHours)<Vco
      Pwind=2.1
  elseif Vci<WindVelocity(1:nHours)<Vr
      Pwind=2.1*((WindVelocity(1:nHours)-Vci)/(Vr-Vci))
  end
  % All Inequality Constraints :-
  powerprob.Constraints.PVout=Ppv<=A*PVeffi*SolarIrradiation'; % Output PV generated power
  powerprob.Constraints.PbCharge=PbCh<=PbChmax*isONbCh; % Battery Limit of allowed Charging 
  powerprob.Constraints.PbDischarge=PbDisch<=PbDischmax*isONbDisch; % Limit of allowed Discharging
  powerprob.Constraints.MaxBatCh=((PbCh(idxHr2Toend,:)*dt)/ec)+(NOMb*SOCb(idxHr2Toend-1,:))<=NOMb; % Maximum battery charge limit
  powerprob.Constraints.NoBChDischtog=isONbCh+isONbDisch<=1 % Forbidden the battery's charging and discharging simultaneously
  powerprob.Constraints.PevCharge=PevCh<=PevChmax*isONevCh; % EV Battery Limit of allowed Charging 
  powerprob.Constraints.PevDischarge=PevDisch<=PevDischmax*isONevDisch; % EV Battery Limit of allowed Discharging
  powerprob.Constraints.MaxEVBatCh=((PevCh(idxHr2Toend,:)*dt)/ec)+(NOMev*SOCev(idxHr2Toend-1,:))<=NOMev; % Maximum EV battery charge limit
  powerprob.Constraints.NoBEVChDischtog=isONevCh+isONevDisch<=1 % Forbidden the ev battery's charging and discharging simultaneously
  % Electric Vehicle's Driving Constraint
  powerprob.Constraints.NOevchout=PevCh(DriveTime)==0; % no ev battery charging when ev is not at home
  powerprob.Constraints.EVdrive=(PevDisch(DriveTime))*dt==2; % Driving electricity demand of ev
  % enforce Battery start=1 when moving from off to on
  % no need to enforce startup=0 at other times since minimizing costs forces it
  powerprob.Constraints.startupConstbCh= -isONbCh(idxHr2Toend-1,:)+isONbCh(idxHr2Toend,:)-startupBch(idxHr2Toend,:)<=0;
  powerprob.Constraints.startupConstbDisch= -isONbDisch(idxHr2Toend-1,:)+isONbDisch(idxHr2Toend,:)-startupBdisch(idxHr2Toend,:)<=0;
  % enforce EV Battery start=1 when moving from off to on
  % no need to enforce startup=0 at other times since minimizing costs forces it
  powerprob.Constraints.startupConstEVch= -isONevCh(idxHr2Toend-1,:)+isONevCh(idxHr2Toend,:)-startupPevCh(idxHr2Toend,:)<=0;
  powerprob.Constraints.startupConstEVdisch= -isONevDisch(idxHr2Toend-1,:)+isONevDisch(idxHr2Toend,:)-startupPevdisCh(idxHr2Toend,:)<=0;
  % % Options for the optimization algorithm, here we set the max time it can run for
  options=optimoptions('intlinprog','Maxtime',1000);
  % Call the optimization solver to find the best solution
  [sol,TotalCost,exitflag,output]=solve(powerprob,'options',options);

The linear constraints and bounds could be inconsistent but I am unable to find the fault. Any help would be greatly appreciated.

Answer 1

[Alan Weiss]

To check whether the bounds and linear constraints are inconsistent, see Converged to an Infeasible Point. Basically, just take your bounds and linear constraints and solve a problem with the objective function 0. For your case, you also have integer constraints, and maybe you could try relaxing them to continuous to start.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Paramvir Singh]

I tried with keeping objective function to zero, relaxing integer variables to continuous and using linprog solver. But result still shows no feasible solution.

[Paramvir Singh]

When I try to run the problem for first interval only instead of 24 intervals, then feasible solution is there. But when I add SOC variable's constraints(State of Charge for batteries) for 24 hour basis then only no feasible solution is shown by the solver. The problem lies with the SOC variable's constraints. Could you please check my code

[ian adrian]

dear Paramvir Singh

I tried to learn about energy management systems at home with the milp method and I got your post. can I see your program and title your paper from IEEE ? for my reference

Answer 2

[awais noor]

% Call the optimization solver to find the best solution

[sol,TotalCost,exitflag,output]=solve(powerprob,'options',options);

Problem faced here an Error as below in matlab:

Error using optim.problemdef.OptimizationProblem/solve

options is not a valid solver. Use 'linprog' or 'intlinprog' instead.

Error in matlbevi (line 115)

[sol,TotalCost,exitflag,output]=solve(powerprob,'options',options);

i m also using this above code for managemnt home i m facing same problem if some one solve it

please give me a solution as soon as possible

([e-mail])

[Alan Weiss]

Please don't ask a question as an answer to a previous question. Instead, start a new question.

As for your question, I suspect that you have an older version of MATLAB, and could find the correct way to pass options in the documentation for your version.

Alan Weiss

MATLAB mathematical toolbox documentation

Answer 3

[ian adrian]

dear Paramvir Singh

I tried to learn about energy management systems at home with the milp method and I got your post. can I see your program and title your paper from IEEE ? for my reference