清空筛选条件
清空筛选条件

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

13 次查看（过去 30 天）
Paramvir Singh
Paramvir Singh 2018-4-20
编辑： ian adrian 2020-7-24
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.

请先登录，再回答此问题。

回答（3 个）

Alan Weiss
Alan Weiss 2018-4-20

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

  3 个评论
显示 1更早的评论
Paramvir Singh
Paramvir Singh 2018-4-21
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
ian adrian 2020-7-24
编辑：ian adrian 2020-7-24
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

请先登录，再进行评论。

awais noor
awais noor 2019-3-13
% 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
(awais2099@live.com)
  1 个评论
Alan Weiss
Alan Weiss 2019-3-13
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

请先登录，再进行评论。

ian adrian
ian adrian 2020-7-24
编辑：ian adrian 2020-7-24
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

类别

在 Help Center 和 File Exchange 中查找有关 Agriculture 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by