MATLAB Answers

Problem-Based Optimization, No feasible solution found

4 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
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
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
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.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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!