Problem with MILP optimizatio problem: "Intlinprog" does not stop
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (1)
Alan Weiss
am 27 Jan. 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
Siehe auch
Kategorien
Mehr zu Get Started with Optimization Toolbox finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!