optimization in smart energy management system.I have used MILP algorithm but i got infeasible solution no point satisfy the constraints.how could i sort out these issues.

14 ビュー (過去 30 日間)
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);
Solving problem using intlinprog. No feasible solution found. Intlinprog stopped because no point satisfies the constraints.

回答 (2 件)

John D'Errico
John D'Errico 2023 年 3 月 27 日
編集済み: John D'Errico 2023 年 3 月 27 日
No solution means exactly what it sounds like. NO solution. And for this class of problem, if it tells you that NO solution exists, it means it. Not just that it could not find a solution, but that none exists. So what do you do?
First, check your equations. CAREFULLY. Look for sign errors. A good thing is that you used an opttimvar approach. This makes your equations much easier to read. So check them all. Then check a second time. Look at the numbers themselves. Did you type in 0.012, when 0.12 was the correct value?
Next, you used rounded numbers here, in essentially ALL of your coefficients. And that is ok. But to me, it suggests we (and by we, I mean YOU) need to look carefully at your equality constraints. The problem is, equality constraints are the most restrictive type of constraint, as they restrict the search space the most. So I would be looking there. The useful trick might be to look at your equality constraints, by replacing them each with a pair of inequality constraints. That is, instead on ONE equality of the form
x + y = 0
replace it with these two constraints:
x + y > -0.1
x + y < 0.1
or some similar appropriate small value. I chose a tolerance of that magnitude, because again, all of your numbers only have 2 significant digits anyway. If no solution is found yet, then expand those tolerances, one at a time. Look to see if suddenly a solution drops out when one of the equalities is loosened significantly.
Then look carefully at that equality. Have you implemented whatever constraint it is properly? Again, is there an error in that constraint?
All of this is not something we can do for you however, as only you know what the various coefficients mean, where they came from. Only you know how you derived the various constraints.
Finally, you should recognize that all of this is just a model. A model represents an approximation of reality. But that does not mean that a solution to everything you write down will exist, and that sometimes the linear approximations in that model might be inappropriate. First though, I would do as I suggest above. And that starts with checking every line you have written there. If necessary, rederive all of those equations, looking for sign errors, or coefficient errors.

Torsten
Torsten 2023 年 3 月 27 日
I have used MILP algorithm but i got infeasible solution no point satisfy the constraints.how could i sort out these issues.
It's hard, but you will have to make a trial and error search within your constraints and see which of them is responsible that your problem has no feasible point.

カテゴリ

Help Center および File ExchangeSources についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by