Problem with MILP optimizatio problem: "Intlinprog" does not stop

11 visualizaciones (últimos 30 días)
lorenzo vallerini
lorenzo vallerini el 27 de En. de 2022
Comentada: Alan Weiss el 27 de En. de 2022
Dear All,
I am currently working on an optimization problem to maximize the revenue from a system consisting of a PV plant. a battery, and an electrolyzer. When I try to solve it, I use "Intlinprog" but the code does not give any resluts and the only thing to do is to force the code to stop ( Ctrl+C). I think the problems is in the constraints, becuase I used "Intlinprog" to solve a different problem and in that case it worked properly.
I also added to "Intlinprog" the option 'MaxTime' but in that case the results present a Relative Gap to high, I can not consder them as the optimal solutions.
Any assistance or insight on how to run "Intlinprog" properly would be greatly appreciated.
nome_simulazione='off_grid';
%% Import PRODUZIONE PV
filename = strcat(pwd , '\PRODUZIONE.csv');
delimiter = {''};
formatSpec = '%f%[^\n\r]';
fileID = fopen(filename,'r');
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'TextType', 'string', 'EmptyValue', NaN, 'ReturnOnError', false);
fclose(fileID);
PRODUZIONE = table(dataArray{1:end-1}, 'VariableNames', {'VarName1'});
clearvars filename delimiter formatSpec fileID dataArray ans;
Ppv = table2array(PRODUZIONE)';
% inserire degradazione annua 0.5%
Ppv1 = [Ppv Ppv*0.995 Ppv*0.99 Ppv*0.985 Ppv*0.980 Ppv*0.975 Ppv*0.97 Ppv*0.965 Ppv*0.960 Ppv*0.955 Ppv*0.95 Ppv*0.945 Ppv*0.94 Ppv*0.935 Ppv*0.93 Ppv*0.925 Ppv*0.920 Ppv*0.915]/1000000; % porto in [MW] da
%% Import PARAMETRI PER DEGRADAZIONE (la tecnologia usata è NMC)
load('aging_param_NMC.mat');
%% TEST VALUES:
PPVp = 100; % potenza nominale di picco del PV [MW]
HHV = 0.0394; % [MWh] potere calorifico di 1 kg di H2
dT = 0.25; % periodo di tempo tra un 'campione' ed il successivo [h]
incremento = 15; % incremento del costo per l'acquisto rispetto al prezzo di vendita
eta_pvg = 0.985*0.995; % rendimento scambio PV-Grid (inverter*traffo)
eta_pvb = 0.985*0.985; % rendimento scambio PV-Battery (inverter*inverter)
eta_be = 0.94; % rendimento scambio Battery-elettrolizzatore (inverter*traffo)
eta_pve=0.94; % rendimento di scambio PV-elettrolizzatore
Pemax=40;
Pbmax=40;
Cnom=8*Pbmax;
Cmin=75;
PriceH2=4;
startCostHIGH =800;
%startCostLOW = 0.2*Pemax;
SOCmin = 10; % valore di SOC minimo consentito
SOCmax = 90; % valore di SOC massimo consentito
SOCini = 50; % SOC iniziale
Price_electricity = 1;
% INITIALIZATIONS:
RicaviGiornalieri = 0;
ore_funz=0;
ore_funz_fine_vita=80000;
tot_day=0;
d = 0;
sum_fd = 0;
sum_fd_cy = 0;
sum_fd_cal = 0;
L = 0;
C_retention = Cnom;
Cnew = Cnom;
a= 0.008439;
b=-0.03224;
c=0.04967;
e=-0.07391;
f=0.997;
eta_battery = a*(Pbmax/Cnom)^4+b*(Pbmax/Cnom)^3+c*(Pbmax/Cnom)^2+e*(Pbmax/Cnom)+f;
eta_c = eta_battery; % rendimento di carica scelto in funzione del Crate massimo
eta_d = eta_c;
%tecnologia NMC (dati per lookup table modello simulink)
RTE = [0.00803226 0.011567 0.015622 0.025834]; % RTE
BP = [0.2 0.3 0.5 1]; % Breakpoints
%% OPTIMIZATION
while ore_funz <= ore_funz_fine_vita && (Cnew/Cnom)*100 >= Cmin
disp('day:')
disp(d+1)
daily_Ppv =(Ppv1(1,(d*96)+1:(d+1)*96)); % [MW] Potenza prodotta dal parco fotovoltaico nel giorno considerato
n = numel(daily_Ppv);
% rendimento elettrolizzatore
yrend = zeros(n,5);
yrend(:,1)=0;
yrend(:,2)=0.66;
yrend(:,3)=0.70;
yrend(:,4)=0.75;
yrend(:,5)=0.72;
% Potenza elettrolizzatore
yPelet=zeros(n,5);
yPelet(:,1)=0.04*Pemax;
yPelet(:,2)=0.25*Pemax;
yPelet(:,3)=0.5*Pemax;
yPelet(:,4)=0.75*Pemax;
yPelet(:,5)=1*Pemax;
% OPTIMIZATION PROBLEM
PowerFlow = optimproblem;
% VARIABLES
Ppvg = optimvar('Ppvg',1,n,'LowerBound',0,'UpperBound',daily_Ppv); % potenza PV -> rete
Ppve = optimvar('Ppve',1,n,'LowerBound',0,'UpperBound',Pemax/eta_pve); % potenza PV -> elettrolizzatore
Ppvb = optimvar('Ppvb',1,n,'LowerBound',0,'UpperBound',Pbmax/eta_pvb);
SOC = optimvar('SOC',1,n,'LowerBound',SOCmin,'UpperBound',SOCmax);
Pbe = optimvar('Pbe',1,n,'LowerBound',0,'UpperBound',Pbmax/eta_d);
Peletlogical = optimvar('Peletlogical',n,11,'Type','integer','LowerBound',0,'UpperBound',1);
z = optimvar('z',n,'Type','integer','LowerBound',0,'UpperBound',1);
Z_X = optimvar('Z_X',1,n,'Type','integer','LowerBound',0,'UpperBound',1); % charging binary variable array
Z_Y = optimvar('Z_Y',1,n,'Type','integer','LowerBound',0,'UpperBound',1); % discharging binary variable array
% OBJECTIVE
RPVG = sum(((Ppvg.*eta_pvg).*Price_electricity)*dT);
RH2 = sum((((Peletlogical(:,1).*yPelet(:,1)).*yrend(:,1))+...
((Peletlogical(:,2).*yPelet(:,2)).*yrend(:,2))+...
((Peletlogical(:,3).*yPelet(:,3)).*yrend(:,3))+...
((Peletlogical(:,4).*yPelet(:,4)).*yrend(:,4))+...
((Peletlogical(:,5).*yPelet(:,5)).*yrend(:,5)))/HHV)*dT*PriceH2;
startingCost = sum((z.*startCostHIGH)*dT);
PowerFlow.Objective = +RPVG-RH2+startingCost;
w = optimexpr(n,1);
idx =2:n;
if d==0
w(1,1)=Peletlogical(1,1);
w(idx,1)=Peletlogical(idx,1)-Peletlogical(idx-1,1)+...
Peletlogical(idx,2)-Peletlogical(idx-1,2)+...
Peletlogical(idx,3)-Peletlogical(idx-1,3)+...
Peletlogical(idx,4)-Peletlogical(idx-1,4)+...
Peletlogical(idx,5)-Peletlogical(idx-1,5);
end
if d>=1
w(1,1)= Peletlogical(1,1)-Results{1,d}.Peletlogical(n,1)+...
Peletlogical(1,2)-Results{1,d}.Peletlogical(n,2)+...
Peletlogical(1,3)-Results{1,d}.Peletlogical(n,3)+...
Peletlogical(1,4)-Results{1,d}.Peletlogical(n,4)+...
Peletlogical(1,5)-Results{1,d}.Peletlogical(n,5);
w(idx,1)=Peletlogical(idx,1)-Peletlogical(idx-1,1)+...
Peletlogical(idx,2)-Peletlogical(idx-1,2)+...
Peletlogical(idx,3)-Peletlogical(idx-1,3)+...
Peletlogical(idx,4)-Peletlogical(idx-1,4)+...
Peletlogical(idx,5)-Peletlogical(idx-1,5);
end
switchcons = w - z <= 0;
% CONSTRAINTS
PowerFlow.Constraints.Ch_on = Ppvb <= (Pbmax/eta_pvb)*Z_X;
PowerFlow.Constraints.Dis_on = Pbe <= (Pbmax/eta_d)*Z_Y;
PowerFlow.Constraints.notPositive = Z_X + Z_Y <= 1;
PowerFlow.Constraints.switchcons = switchcons;
PowerFlow.Constraints.tot=Peletlogical(:,1)+...
Peletlogical(:,2)+...
Peletlogical(:,3)+...
Peletlogical(:,4)+...
Peletlogical(:,5)<=1;
PowerFlow.Constraints.tot2=((Peletlogical(:,1).*yPelet(:,1)))+...
((Peletlogical(:,2).*yPelet(:,2)))+...
((Peletlogical(:,3).*yPelet(:,3)))+...
((Peletlogical(:,4).*yPelet(:,4)))+...
((Peletlogical(:,5).*yPelet(:,5)))-Pbe'== +Ppve';
PowerFlow.Constraints.PPV = daily_Ppv == Ppvb + Ppvg + Ppve;
PowerFlow.Constraints.rampa=z+...
Peletlogical(:,1)*2+...
Peletlogical(:,2)*2+...
Peletlogical(:,3)+...
Peletlogical(:,4)*2+...
Peletlogical(:,5)*2<=2;
PowerFlow.Constraints.SOC = optimconstr(n); % vincoli sullo Stato di Carica
for j=1:n
if j==1
PowerFlow.Constraints.SOC = SOC(j) == SOCini;
else
PowerFlow.Constraints.SOC(j) = SOC(j) == SOC(j-1)+((Ppvb(j)*eta_pvb)-(Pbe(j)*eta_d))/Cnew*(100*dT);
end
end
PowerFlow.Constraints.SOC2= SOC(end)==50;
options = optimoptions('intlinprog','Display','final','MaxTime',7);
[Result,fval,exitflag,output]= solve(PowerFlow,'options',options);
% RESULTS SAVING
SOCini = Result.SOC(length(SOC));
Results(d+1) = {Result}; % salvataggio dei profili di potenza giornalieri ottenuti dall'ott.
RicaviGiornalieri(d+1) = -fval;
index = zeros(1,length(Result.Peletlogical));
index(Result.Peletlogical >= 0.1)=1;
ore_funz = ore_funz+sum(index)/4;
tot_day=tot_day+length(daily_Ppv)/96;
years_of_life = (tot_day)/365;
d=d+1;
% BATTERY DEGRADATION:
% BATTERY TEMPERATURE EVALUATION (Simulink)
t = 0:0.25:23.75; % vettore dei tempi per il giorno
P_batt = timeseries((Result.Pbe-Result.Ppvb),t*3600); % timeseries per Simulink [MW]
T_amb = 22; % [°C]
%dati per tecnologia NMC
n_module = Cnom/0.013; % numero di moduli utilizzati
s = (1.162*0.110*2)+(0.445*0.110*2); % superficie di scambio per ogni modulo [m^2]
cell_area = n_module*s; % [m^2]
h_conv = 200; % heat transfer coefficient (convezione forzata)
w = 89; % [kg] peso cella
cell_mass = n_module*w;
cell_Cp_heat = 1600; % specific heat [J/kg/K]
tic
t_simulation = t(length(t))*3600;
if d==0
T_init = T_amb+273.15;
end
if d>=1
T_init = T_bat.Data(end)+273.15; % in ogni giorno T ini. = T fin. giorno precedente
end
sim('modello_Temperatura_Batteria')
toc
% fattori di stress:
% tecnologia NMC
[fd_aging_cy,fd_aging_cal,S_temp,Temp] = battery_degradation(Result,aging_param_NMC,T_bat,t);
sum_fd_cy = sum_fd_cy+fd_aging_cy;
sum_fd_cal = sum_fd_cal+fd_aging_cal;
sum_fd = sum_fd+fd_aging_cy+fd_aging_cal;
L(d+1) = 1-aging_param_NMC.alfasei*exp(-aging_param_NMC.betasei*sum_fd)-(1-aging_param_NMC.alfasei)*exp(-sum_fd);
Cnew = Cnom*(1-L(d+1)); % nuovo valore di capacità per il giorno successivo
C_retention(1,d+2) = Cnew; % profilo di variazione della capacità
disp('nuovo valore di capacità:')
disp(Cnew)
end
save('off_grid_C/D')
  4 comentarios
Torsten
Torsten el 27 de En. de 2022
So the solver does not find a feasible point to start with ?
lorenzo vallerini
lorenzo vallerini el 27 de En. de 2022
Editada: lorenzo vallerini el 27 de En. de 2022
I don't know if this could help, but this is what I get when i run the code(with Branch and Bound it continue to explore nodes, I stoped the program after 20 minutes and it was sill going):
Solving problem using intlinprog.
LP: Optimal objective value is -35440.200417.
Cut Generation: Applied 8 Gomory cuts,
35 implication cuts, 17 clique cuts,
11 flow cover cuts, 16 mir cuts,
and 3 cover cuts.
Lower bound is -35155.738375.
Heuristics: Found 1 solution using total rounding.
Upper bound is -34964.284331.
Relative gap is 0.55%.
Cut Generation: Applied 2 Gomory cuts, and 2 clique cuts.
Lower bound is -35155.738365.
Relative gap is 0.55%.
Branch and Bound:
nodes total num int integer relative
explored time (s) solution fval gap (%)
2704 4.99 2 -3.507074e+04 2.423465e-01
12704 14.22 2 -3.507074e+04 2.423465e-01
22704 22.71 2 -3.507074e+04 2.423465e-01

Iniciar sesión para comentar.

Respuestas (1)

Alan Weiss
Alan Weiss el 27 de En. de 2022
You might notice that the default value of the MaxTime option is 7200 seconds. That is because it is typical for intlinprog to take hours to solve a problem. Sorry, that is just the way it is for integer linear programming solvers.
For advice on how to speed intlinprog, see Tuning Integer Linear Programming. But don't have the wrong expectations: no matter what you do, solving any particular problem can take a very long time.
Alan Weiss
MATLAB mathematical toolbox documentation
  5 comentarios
lorenzo vallerini
lorenzo vallerini el 27 de En. de 2022
Is it possible to speed up the optimization by changing some of the constraints I put? Because in a different optimization problem (where instead of the battery there was the electric grid) intlinprog in a matter of seconds, finding the optimal solution (IntegerTolerance = 1e-05,the default value).
Alan Weiss
Alan Weiss el 27 de En. de 2022
I repeat: For advice on how to speed intlinprog, see Tuning Integer Linear Programming.
Alan Weiss
MATLAB mathematical toolbox documentation

Iniciar sesión para comentar.

Categorías

Más información sobre Get Started with Optimization Toolbox en Help Center y File Exchange.

Productos


Versión

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by