Problem-Based Optimization, No feasible solution found

5 views (last 30 days)
In this simplified version of a system designed to produce hydrogen from pv, I want to find out the lowest possible capacity of the electrolyzer?
The problem is that when I run the model it says that there is no solution avaliable. However, when I fill in the capacity of the electrolyzer by hand I can easily find a solution for the constraints.
%% Problem-Based Optimization
load('CF_pv.mat'); % Hourly data of capacity factor sun
energyprob = optimproblem;
PV = 748; % in MW
Electrolyzer = optimvar('Electrolyzer',1,'LowerBound',0); % Capacity of ELectrolyzer in MW
Electrolyzer_costs =[(200+0.7*200)];
costs = Electrolyzer_costs*Electrolyzer;
%Objective Function
energyprob.Objective = costs;
%% Optimization constraints
PV_Load = PV*CF_pv;
T = 8760; % hours in a year
Bat_In = optimvar('Bat_In',T,'LowerBound',0); % Electricity send to battery
Bat_Out = optimvar('Bat_Out',T,'LowerBound',0); % Electricity send from battery to Electrolyzer
Bat_In_c = optimconstr(T);
Bat_Out_c = optimconstr(T);
for t=1:T
Bat_In_c(t) = Bat_In(t) == (PV_Load(t)-1.1*Electrolyzer); % Abundant Electricity, above 110% of the limit of the electrolyzer
Bat_Out_c(t) = Bat_Out(t) == 0.1*Electrolyzer-PV_Load(t); % Electricity necessary to keep electrolyzer running
end
Elec_In = PV_Load - Bat_In + Bat_Out; % Generated electricity - input in battery + electricty from battery to electrolyzer
energyprob.Constraints.Bat_In_c = Bat_In_c;
energyprob.Constraints.Bat_Out_c = Bat_Out_c;
Desired_Elec_In = 1.7618*10^6; % Desired electricty input for a whole year
cons1 = Elec_In == sum(Desired_Elec_In);
energyprob.Constraints.cons1 = cons1;
[sol,fval] = solve(energyprob);
sol.Electrolyzer % display electrolyzer capacity

Accepted Answer

Alan Weiss
Alan Weiss on 2 Jun 2021
Your problem is infeasible as given. Using the routines in Investigate Linear Infeasibilities, I found the following. Convert the problem to LP matrices:
problem = prob2struct(energyprob);
Add a null Aineq and bineq inequality constraint.
Aineq = zeros(1,size(problem.Aeq,2));
bineq = 0; % The code needs an Aineq argument of the correct size
I then ran the code in Remove iIS In A Loop, and found that the following are conflicting lines in the problem.Aeq matrix together with the problem.beq data (conflict with the bounds):
Row 8760 and 17513
Row 8759 and 17512
Row 8758 and 17511
etc.
It just kept spitting them out, so I stopped the process. Armed with this knowledge I am sure that you can find the error in your setup.
By the way, I recently found some errors in the published code for the generatecover function. The corrections you must make in the generatecover code:
% Look for this line
[ineq_iis,eq_iis,ncalls,fval] = elasticfilter(Aineq(activeA,:),bineq(activeA),Aeq(activeAeq),beq(activeAeq),lb,ub);
% change it to
[ineq_iis,eq_iis,ncalls,fval] = elasticfilter(Aineq(activeA,:),bineq(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub);
% Look for this line
[ineq_iis,eq_iis,ncalls,fval] = elasticfilter(Aineq(activeA),bineq(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub);
% change it to
[ineq_iis,eq_iis,ncalls,fval] = elasticfilter(Aineq(activeA,:),bineq(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub);
But I suggest that you do not run generatecover on this problem; I have been running it for hours and it is still not finished.
Good luck,
Alan Weiss
MATLAB mathematical toolbox documentation
  2 Comments
Alan Weiss
Alan Weiss on 2 Jun 2021
Just a bit more info. Looking at Rows 8760 and 17513, each has only two nonzero entries. 8760:
Aeq(8760,8760) = 1.0000
Aeq(8760,17521) = 1.1000
beq(8760) = 0
Together, these mean that x(8760) + 1.1*x(17521) = 0. But the lower bound on all x components is 0, so this implies x(8760) = x(17521) = 0.
Examine row 17513.
Aeq(17513,17513) = 1
Aeq(17513,17521) = -0.1
beq(17513) = -115.8937
These mean x(17513) - 0.1*x(17521) = -115.8937.
Therefore we must have x(17521) > 0 (actually x(17521) > 1158.937). But we already found that x(17521) = 0. So this is a contradiction.
Alan Weiss
MATLAB mathematical toolbox documentation
Rogier Doodeman
Rogier Doodeman on 9 Jun 2021
Thank you for your answer, this was indeed the problem.
For those interested the code delvired the right answer when I removed the contraints determining the Bat_In and Bat_Out value and instead wrote contraints 2 and 3.
%% Problem-Based Optimization
load('CF_pv.mat'); % Hourly data of capacity factor sun
energyprob = optimproblem('ObjectiveSense','minimize');
PV = optimvar('PV','Type','integer','LowerBound',500); % Capacity of pv panels in MW
Electrolyzer = optimvar('Electrolyzer','Type','integer','LowerBound',0); % Capacity of ELectrolyzer in MW
PV_costs =[(550+5.5)*10^3]*PV; %[(550+5.5)*(max(Output_Solar_DC)] times the peak output so the price is per meter
Electrolyzer_costs =(200+0.7*200)*Electrolyzer;
%% Optimization constraints
PV_Load = PV*CF_pv;
Desired_Elec_In = 1.7618*10^6; % Desired electricty input for a whole year
T = 8760;
M = 500000; % Value that can not be surpassed
Bat_In = optimvar('Bat_In',T,'LowerBound',0,'UpperBound',M);
Bat_Out = optimvar('Bat_Out',T,'LowerBound',-M,'UpperBound',0);
Z_X = optimvar('Z_X',T,'Type','integer','LowerBound',0,'UpperBound',1); % charging binary variable array
Z_Y = optimvar('Z_Y',T,'Type','integer','LowerBound',0,'UpperBound',1); % discharging binary variable array
Bat_Mut = Bat_In + Bat_Out; % Battery mutation per hour, Bat_Out are all negative values so it will subtract
Elec_In = PV_Load - Bat_Mut;
%Objective Function
energyprob.Objective = PV_costs + Electrolyzer_costs +sum(Z_X) + sum(Z_Y); % + Battery_Ene_costs + Battery_Pow_costs
cons1 = sum(Elec_In) == Desired_Elec_In;
cons2 = Elec_In <= 1.1*Electrolyzer;
cons3 = Elec_In >= 0.1*Electrolyzer;
cons4 = Bat_In <= M*Z_X;
cons5 = Bat_Out >= -M*Z_Y;
cons6 = Z_X + Z_Y <= 1;
energyprob.Constraints.cons1 = cons1;
energyprob.Constraints.cons2 = cons2;
energyprob.Constraints.cons3 = cons3;
energyprob.Constraints.cons4 = cons4;
energyprob.Constraints.cons5 = cons5;
energyprob.Constraints.cons6 = cons6;
[sol,fval] = solve(energyprob)
sol.PV
sol.Electrolyzer % display electrolyzer capacity
To make sure the battery does not charge and discharge at the same time constraints 4 - 6 are used. For more information on this you can look at this thread. https://nl.mathworks.com/matlabcentral/answers/577642-energy-storage-optimisation-problem-separate-charge-and-discharge-constraint

Sign in to comment.

More Answers (0)

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!