admin管理员组文章数量:1244316
I'm trying to solve a MILP problem in Matlab Optimization Toolbox. I'm encountering with problems with the product of two decision variable of my problem which are bounded.
- 0<= x <= 60
- 0 <= y <= 8760
The product is inside the objective function that is minimization of the cost. For this reason when I tried with the method .pdf explained here, I obtained a value always equal to zero, that is not what I want. I'd like to a value in time that is equal to the product x*y. The variable x is a non linear function of time which represents a power, while y is the cumulative sum of the working hours of operation of the system, that goes from 0, if the system is not working, to 8760 if the system work every hour.
clear all
close all
clc
addpath('C:\gurobi1103\win64\examples\matlab')
%% Input data
% Define the main path manually, e.g. 'C:\Users\My_models'
path_main = 'C:\Users\nanf\Desktop\H2-districts';
% path for input files
path_input = fullfile(path_main, 'Input');
% Define the path for functions
path_functions = fullfile(path_main, 'Functions');
% Reading an Excel file from the Input directory
excelFilePath1 = fullfile(path_input, 'Inputs_BoundaryLoads_2022.xlsx');
excelFilePath2 = fullfile(path_input, 'nestel_demand.csv');
excelFilePath3 = fullfile(path_input, 'PriceImpExp_ele.xlsx');
Input1 = readtable(excelFilePath1);
Input2 = readtable(excelFilePath2);
Input3 = readtable(excelFilePath3);
%% Pre-processing
nHours = 8760; % Number of hours simulated [8760h]
irradiance = Input1.G; % Hourly solar irradiance [W/m^2]
idxHr2ToEnd = (2:nHours)'; % Hours until the end
Time = (1:nHours)'; % Time vector
days = nHours/24; % Number of days simulated
weeks = days/7; % Number of weeks simulated
clear t;
linew = 1.5;
font = 18;
%% INPUT PARAMETERS
% general inputs
bigM = 1e8; % Large number for big-M constraints
MaxSimTime = 900; % Maximum time for MILP solver [s]
deltat = 3600; % time step di 1h [s]
% Electrical district demand
P_load = Input2.P_el_nest_h; % NEST electrical demand [kW]
% Fixed technical parameters
HHV = 39.39 * 3.6 * 10^3; % Hydrogen higher heating value [kJ/kg]
% Costs and Efficencies
c_h2 = 2; % cost of hydrogen [CHF/kg]
c_gridimp = Input3.Priceimpgreen; % [CHF/kWh]
c_gridexp = Input3.PriceexpPV; % [CHF/kWh]
eff_b = 0.95; % constant efficiency of the battery (round trip efficency)
eff_PV = 0.17; % constant efficiency of the PV system
eff_FC = 0.50; % Nominal electrical efficency of FC
% Lifetime in hours for PEMFC from Technology RoadMap-Hyydrogen and Fuel Cells,IEA(2015)
nH_lifetime = 40000; % [h]
% PEM fuel cell lifetime was defined by the time elapsed until 10% of the initial voltage/performance is lost
Delta_eff_FC = 0.1*eff_FC; % [-]
% Fuel Cell degradation rate as a constant value in Nani-Model
Const1 = Delta_eff_FC/nH_lifetime; % [1/h]
% further constant for MILP model from Taylor expansion
Const3 = Const1/HHV/eff_FC^2; %[kg/h/kJ]
% Const2 = Const3*P_FC_opt;
% Fuel Cell (FC)
S_FC_min = 0; % Minimum size FC [kW]
S_FC_max = 60; % Maximum size FC [kW]
% PV
P_PV_peak_ref = 100; % Reference value peak power PV [kW]
Area_PV_ref = P_PV_peak_ref/eff_PV; % Reference area PV [m2]
Area_PV_min = 0; % Minimum area PV [m2]
Area_PV_max = Area_PV_ref*2; % Maximum area PV [m2]
% Battery
C_b_max = 96 * 3.6 * 10^3; % Battery capacity [kJ]
C_b_min = C_b_max / 3; % MIN Battery capacity [kJ]
P_b_max = C_b_max / 3600; % Battery capacity (1C rate) [kWh]
% unit prices and lifetime components
d = 0.05;%0.04; % Discount rate
ann = d / (1 - (1 + d)^(-20)); % annuity factor calculated with plant lifetime
UP_PV = 1300;%1240; % Unit price PV [CHF/kW_p]
life_PV = 25;%30; % Lifetime PV [years]
maint_PV = 0.01;%0.0158; % Annual cost maintenance PV, frac total cost
ann_PV = d / (1 - (1 + d)^(-life_PV));
UP_b = 1000;%600 % Unit price battery [CHF/kWh]
life_b = 12;%10; % Lifetime battery [years]
maint_b = 0.02; % frac capex
ann_b = d / (1 - (1 + d)^(-life_b));
UP_FC = 950;%3000 % Unit price FC [CHF/kW]
life_FC = 7.5;%15 % Lifetime FC [years]
maint_FC = 0.024;%0.02 % Annual cost maintenance FC, frac capex cost
ann_FC = d / (1 - (1 + d)^(-life_FC));
aux_maxs = (S_FC_max+nHours);
aux_mins = 0;
aux_maxd = S_FC_max;
aux_mind = -nHours;
%% PWL
run PWL_productcontvariables.m
%% Define the optimization problem and the optimization variables
sizingprob = optimproblem;
%% DECISION VARIABLES
% SIZING decision variables
% Battery capacity in [kJ]
C_b = optimvar('C_b','LowerBound',C_b_min,'UpperBound',C_b_max); % C_b actual installed capacity,C_b_max= max capacity we can have
% Fuel Cell size in [kW]
S_FC = optimvar('S_FC','LowerBound',S_FC_min,'UpperBound',S_FC_max);
% PV area in [m2]
Area_PV = optimvar('Area_PV','LowerBound',Area_PV_min,'UpperBound',Area_PV_max);
% OPERATIONAL decision variables
% fuel cell generated power in [kW]
P_FC = optimvar('P_FC',nHours,'LowerBound',0,'UpperBound',S_FC_max);
% additional variables for FC operation implementation
P_FC_On = optimvar('P_FC_On',nHours,'LowerBound',0,'UpperBound',S_FC_max);
delta_On = optimvar('FC_On',nHours,'Type','integer','LowerBound',0,'UpperBound',1);
% Degradation auxiliary variable [h]
w_hours = optimvar('w_hours',nHours,'LowerBound',0,'UpperBound',nHours); %cumulative summation of FC_On
% mass flow rate of H2 [kg/h]
m_flow_H2 = optimvar('m_flow_H2',nHours,'LowerBound',0,'UpperBound',S_FC_max/eff_FC/HHV.*3600*3);
% imported power from the grid in [kW]
P_imp = optimvar('P_imp',nHours,'LowerBound',0);
% exported power to the grid in [kW]
P_exp = optimvar('P_exp',nHours,'LowerBound',0);
% Battery energy content in [kWh]
E_b = optimvar('E_b',nHours,'LowerBound',0,'UpperBound',C_b_max);
% Battery charging power in [kW]
P_b_ch = optimvar('P_b_ch',nHours,'LowerBound',0,'UpperBound',P_b_max);
% Battery discharging power in [kW]
P_b_disch = optimvar('P_b_disch',nHours,'LowerBound',0,'UpperBound',P_b_max);
%% Additional Decision Variables for PWA Approximation
% Squared variables for X and Y
square_X = optimvar('square_X', nHours,'LowerBound', 0,'UpperBound',(nHours+S_FC_max)^2*2); % f(t)=x^2, x=(P_fc+w_hours)^2
square_Y = optimvar('square_Y', nHours, 'LowerBound', 0, 'UpperBound',(nHours+S_FC_max)^2*2); % f(t)=y^2, y=(P_fc-w_hours)^2
X_var = optimvar ('X_var',nHours, 'LowerBound',aux_mins,'UpperBound',aux_maxs);
Y_var = optimvar ('Y_var',nHours,'LowerBound',aux_mind,'UpperBound',aux_maxd);
delta_sum1 = optimvar('delta_sum1','Type','integer','LowerBound',0,'UpperBound',1);
delta_sum2 = optimvar('delta_sum2','Type','integer','LowerBound',0,'UpperBound',1);
delta_sum3 = optimvar('delta_sum3','Type','integer','LowerBound',0,'UpperBound',1);
delta_sum4 = optimvar('delta_sum4','Type','integer','LowerBound',0,'UpperBound',1);
aux_sum1 = optimvar('aux_sum1','LowerBound',aux_mins,'UpperBound',aux_maxs);
aux_sum2 = optimvar('aux_sum2','LowerBound',aux_mins,'UpperBound',aux_maxs);
aux_sum3 = optimvar('aux_sum3','LowerBound',aux_mins,'UpperBound',aux_maxs);
aux_sum4 = optimvar('aux_sum4','LowerBound',aux_mins,'UpperBound',aux_maxs);
delta_diff1 = optimvar('delta_diff1','Type','integer','LowerBound',0,'UpperBound',1);
delta_diff2 = optimvar('delta_diff2','Type','integer','LowerBound',0,'UpperBound',1);
delta_diff3 = optimvar('delta_diff3','Type','integer','LowerBound',0,'UpperBound',1);
delta_diff4 = optimvar('delta_diff4','Type','integer','LowerBound',0,'UpperBound',1);
aux_diff1 = optimvar('aux_diff1','LowerBound',aux_mind,'UpperBound',aux_maxd);
aux_diff2 = optimvar('aux_diff2','LowerBound',aux_mind,'UpperBound',aux_maxd);
aux_diff3 = optimvar('aux_diff3','LowerBound',aux_mind,'UpperBound',aux_maxd);
aux_diff4 = optimvar('aux_diff4','LowerBound',aux_mind,'UpperBound',aux_maxd);
% Reformulated product term
PFC_w_hours = optimvar('PFC_w_hours', nHours, 'LowerBound', 0, 'UpperBound', S_FC_max*nHours);
%% Derived variables
% PV generated power
P_PV = irradiance.*eff_PV*Area_PV/1000; % [kW]
P_PV_peak = 1000*eff_PV*Area_PV/1000; % [kW]
% C-rate of battery [kWh]
P_b_lim = C_b / 3600; % new max capacity in [kWh]
%% CONSTRAINTS
% energy balance
sizingprob.Constraints.EnBalance = (P_FC + P_PV + P_b_disch + P_imp) == (P_b_ch + P_load + P_exp);
% BATTERY
% sizingprob.Constraints.NoSimultaneousChDisch = discharging_on + charging_on <= ones(nHours,1);
sizingprob.Constraints.PowerBatt_ch_0 = P_b_ch(1) == 0;
sizingprob.Constraints.PowerBatt_disch_0 = P_b_disch(1) == 0;
%sizingprob.Constraints.Ch_on1 = P_b_ch <= P_b_max * charging_on;
sizingprob.Constraints.Ch_on2 = P_b_ch <= P_b_lim;
%sizingprob.Constraints.Disch_on1 = P_b_disch <= P_b_max * discharging_on;
sizingprob.Constraints.Disch_on2 = P_b_disch <= P_b_lim;
sizingprob.Constraints.E_b_update = E_b(idxHr2ToEnd) - E_b(idxHr2ToEnd-1) == P_b_ch(idxHr2ToEnd)*deltat - P_b_disch(idxHr2ToEnd)*deltat;
sizingprob.Constraints.E_b_cont = E_b(1) == E_b(end);
sizingprob.Constraints.E_b_max = E_b <= C_b;
% sizingprob.Constraints.E_b_min = E_b >= 0.2 * C_b;
% FC
sizingprob.Constraints.MaxPowerFC = P_FC <= P_FC_On;
% sizingprob.Constraints.MinPowerFC = P_FC >= 0.2 * P_FC_On;
sizingprob.Constraints.PowerFC1 = P_FC_On <= S_FC_max * delta_On;
sizingprob.Constraints.PowerFC2 = P_FC_On <= S_FC;
sizingprob.Constraints.PowerFC3 = P_FC_On >= S_FC - S_FC_max.*(ones(nHours,1) - delta_On);
sizingprob.Constraints.PowerFC4 = P_FC <= S_FC_max * delta_On;
% Constraints for cumulative summation
sizingprob.Constraints.SumCum_0 = w_hours(1) == delta_On(1);
sizingprob.Constraints.SumCum_update = w_hours(idxHr2ToEnd) == w_hours(idxHr2ToEnd-1) + delta_On(idxHr2ToEnd);
sizingprob.Constraints.pwl_segmentactsum = delta_sum1+delta_sum2+delta_sum3+delta_sum4==1 ;
sizingprob.Constraints.pwl_segmentactdiff = delta_diff1+delta_diff2+delta_diff3+delta_diff4==1 ;
%% (Pfc+whours)^2
% for auxiliary variable
sizingprob.Constraints.auxsum11 = aux_sum1 <= aux_maxs*delta_sum1;
sizingprob.Constraints.auxsum12 = aux_sum1 >= aux_mins*delta_sum1;
sizingprob.Constraints.auxsum13 = aux_sum1 <= X_var;
sizingprob.Constraints.auxsum14 = aux_sum1 >= X_var - aux_maxs *(1-delta_sum1);
sizingprob.Constraints.auxsum21 = aux_sum2 <= aux_maxs*delta_sum2;
sizingprob.Constraints.auxsum22 = aux_sum2 >= aux_mins*delta_sum2;
sizingprob.Constraints.auxsum23 = aux_sum2 <= X_var;
sizingprob.Constraints.auxsum24 = aux_sum2 >= X_var - aux_maxs *(1-delta_sum2);
sizingprob.Constraints.auxsum31 = aux_sum3 <= aux_maxs*delta_sum3;
sizingprob.Constraints.auxsum32 = aux_sum3 >= aux_mins*delta_sum3;
sizingprob.Constraints.auxsum33 = aux_sum3 <= X_var;
sizingprob.Constraints.auxsum34 = aux_sum3 >= X_var - aux_maxs *(1-delta_sum3);
sizingprob.Constraints.auxsum41 = aux_sum4 <= aux_maxs*delta_sum4;
sizingprob.Constraints.auxsum42 = aux_sum4 >= aux_mins*delta_sum4;
sizingprob.Constraints.auxsum43 = aux_sum4 <= X_var;
sizingprob.Constraints.auxsum44 = aux_sum4 >= X_var - aux_maxs *(1-delta_sum4);
% sizingprob.Constraints.z_plus_lb_seg1 = X_var >= xxs(1) * delta_sum1;
sizingprob.Constraints.z_plus_ub_seg1 = X_var <= xxs(2) * delta_sum1 + xxs(3) * delta_sum2 + xxs(4)*delta_sum3+xxs(5)*delta_sum4;
% sizingprob.Constraints.z_plus_lb_seg2 = X_var >= xxs(2) * delta_sum2;
% sizingprob.Constraints.z_plus_ub_seg2 = X_var <= xxs(3) * delta_sum2;
% sizingprob.Constraints.z_plus_lb_seg3 = X_var >= xx(3) * delta_sum3;
% sizingprob.Constraints.z_plus_ub_seg3 = X_var <= xx(4) * delta_sum3;
% sizingprob.Constraints.z_plus_lb_seg4 = X_var >= xx(4) * delta_sum4;
% sizingprob.Constraints.z_plus_ub_seg4 = X_var <= xx(5) * delta_sum4;
sizingprob.Constraints.V_PWA_X = square_X == aas(1)* aux_sum1 + bbs(1) * delta_sum1 + aas(2) * aux_sum2 + bbs(2)* delta_sum2 +aas(3) * aux_sum3 + bbs(3)* delta_sum3 + aas(4) * aux_sum4 + bbs(4)* delta_sum4;
%% DIFF
% for auxiliary variable
sizingprob.Constraints.auxdiff11 = aux_diff1 <= aux_maxd*delta_diff1;
sizingprob.Constraints.auxdiff12 = aux_diff1 >= aux_mind*delta_diff1;
sizingprob.Constraints.auxdiff13 = aux_diff1 <= Y_var;
sizingprob.Constraints.auxdiff14 = aux_diff1 >= Y_var - aux_maxd *(1-delta_diff1);
sizingprob.Constraints.auxdiff21 = aux_diff2 <= aux_maxd*delta_diff2;
sizingprob.Constraints.auxdiff22 = aux_diff2 >= aux_mind*delta_diff2;
sizingprob.Constraints.auxdiff23 = aux_diff2 <= Y_var;
sizingprob.Constraints.auxdiff24 = aux_diff2 >= Y_var - aux_maxd *(1-delta_diff2);
sizingprob.Constraints.auxdiff31 = aux_diff3 <= aux_maxd*delta_diff3;
sizingprob.Constraints.auxdiff32 = aux_diff3 >= aux_mind*delta_diff3;
sizingprob.Constraints.auxdiff33 = aux_diff3 <= Y_var;
sizingprob.Constraints.auxdiff34 = aux_diff3 >= Y_var - aux_maxd *(1-delta_diff3);
sizingprob.Constraints.auxdiff41 = aux_diff4 <= aux_maxd*delta_diff4;
sizingprob.Constraints.auxdiff42 = aux_diff4 >= aux_mind*delta_diff4;
sizingprob.Constraints.auxdiff43 = aux_diff4 <= Y_var;
sizingprob.Constraints.auxdiff44 = aux_diff4 >= Y_var - aux_maxd *(1-delta_diff4);
% sizingprob.Constraints.z_minus_lb_seg1 = Y_var >= xxd(1) * delta_diff1;
sizingprob.Constraints.z_minus_ub_seg1 = Y_var <= xxd(2) * delta_diff1+xxd(3) * delta_diff2+xxd(4)*delta_diff3+xxd(5)*delta_diff4;
% sizingprob.Constraints.z_minus_lb_seg2 = Y_var >= xxd(2) * delta_diff2;
% sizingprob.Constraints.z_minus_ub_seg2 = Y_var <= xxd(3) * delta_diff2;
% sizingprob.Constraints.z_minus_lb_seg3 = Y_var >= xx(3) * delta_diff3;
% sizingprob.Constraints.z_minus_ub_seg3 = Y_var <= xx(4) * delta_diff3;
% sizingprob.Constraints.z_minus_lb_seg4 = Y_var >= xx(4) * delta_diff4;
% sizingprob.Constraints.z_minus_ub_seg4 = Y_var <= xx(5) * delta_diff4;
sizingprob.Constraints.V_PWA_Y = square_Y == aad(1) * aux_diff1 + bbd(1) * delta_diff1 + aad(2) * aux_diff2 + bbd(2) * delta_diff2 + aad(3) * aux_diff3 + bbd(3) * delta_diff3 + aad(4) * aux_diff4 + bbd(4)* delta_diff4;
% Reformulated product constraint using both PWA functions
sizingprob.Constraints.ProductApprox1 = PFC_w_hours == (1/4) * (square_X - square_Y);
sizingprob.Constraints.ProductApprox2 = PFC_w_hours <=nHours*S_FC_max;
sizingprob.Constraints.ProductApprox3 = PFC_w_hours >=0;
% mass flow rate calculation
sizingprob.Constraints.massflowH2= m_flow_H2 == (P_FC./eff_FC./HHV).*3600 + PFC_w_hours*Const3*3600; % [kg/h]
%% OBJECTIVE FUNCTION
cost_inst = (P_PV_peak * UP_PV * ann_PV + S_FC * UP_FC * ann_FC + C_b/3600 * UP_b * ann_b)/1000; % [kEUR/y]
cost_imp = sum(c_h2 * m_flow_H2) / 1000 + sum(P_imp .* c_gridimp) / 1000; % [kEUR/y]
cost_exp = sum(P_exp .* c_gridexp) / 1000; % [kEUR/y]
cost_maint = (maint_PV * P_PV_peak * UP_PV + maint_FC * S_FC * UP_FC + maint_b * UP_b * C_b/3600)/1000; % [kEUR/y]
cost = cost_inst + cost_imp - cost_exp + cost_maint;
% set objective
sizingprob.Objective = cost;
%% Solve optimization problem
intcon = [];
[solution,fval,reasonSolverStopped] = solve(sizingprob);
本文标签： matlabHow to linearize the product of two continuous variablesStack Overflow
版权声明：本文标题：matlab - How to linearize the product of two continuous variables? - Stack Overflow 内容由网友自发贡献，该文观点仅代表作者本人， 转载请联系作者并注明出处：http://www.betaflare.com/web/1740203296a2240600.html， 本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容，一经查实，本站将立刻删除。
发表评论