Problem-Based Optimization, No feasible solution found
5 views (last 30 days)
Show older comments
Rogier Doodeman
on 31 May 2021
Commented: Rogier Doodeman
on 9 Jun 2021
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
0 Comments
Accepted Answer
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
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
More Answers (0)
See Also
Categories
Find more on Nonlinear Optimization in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!