How can I solve a large-size MILP problem?
Hello everyone,
I try to solve a large-size MILP problem as the figure shows:
In terms of solving, I think it faces two main challenges:
(1) This is a NP-hard problem due to some w are binary. Besides, it somehow like a sample-based stochastic problem, which easily leads to the curse of dimension.
(2) The variable B couple every sub-MILP problem, and its norm-1 term may result some trouble.
With n = 274 (the model involves 274 scenarios), I carry out some testing. Gurobi reports me that it requires several hours for root relaxation, so I think it is too difficult to solve if I do not prepocess the model (even though I have done a warm-starting and use a PC with i9 8-cores CPU and 32GB RAM).
In terms of prepocessing, could you please share some ideas with me? How about Lagrange-relaxation-based decomposition to the coupled constraints Bx? Or any other advice for tuning the Gurobi parameters?
My code is as following, it is based on MATLAB-YALMIP-Gurobi to implement. Any advice to improve the code is highly appreciated.
DO_UC_test_database;
Number_scenarios_train = 274;
load Recorder_obj_train_test;
load Realization_load_bus_train_test;
load Feature_D_train_test;
load Recorder_S_onoff;
load Recorder_S_su;
load Recorder_S_sd;
z = Recorder_obj_train;
lamda_set = [10000, 100000, 1000000, 1000000];
lamda_no = 2;
lamda = lamda_set(lamda_no);
B = sdpvar(Number_horizon*Number_RES,Number_horizon+Number_horizon*Number_RES);
Constraint = [];
% t for t
% Section #1
Constraint = Constraint + [B(001:002,002:024) == 0] + [B(001:002,027:072) == 0];
for i = 2:23
Constraint = Constraint + [B(Number_RES*(i-1)+1:Number_RES*i,001:i-1) == 0] + [B(Number_RES*(i-1)+1:Number_RES*i,i+1:024) == 0]...
+ [B(Number_RES*(i-1)+1:Number_RES*i,025:024+Number_RES*(i-1)) == 0] + [B(Number_RES*(i-1)+1:Number_RES*i,024+i*Number_RES+1:072) == 0];
end
% Section #24
Constraint = Constraint + [B(047:048,001:023) == 0] + [B(047:048,025:070) == 0];
%% ---------------------------binary variabels----------------------------%
Decision_State_onoff = cell(Number_scenarios_train,1);
Decision_State_startup = cell(Number_scenarios_train,1);
Decision_State_shutdown = cell(Number_scenarios_train,1);
for NST = 1:Number_scenarios_train
Decision_State_onoff{NST} = binvar(Number_gen,Number_horizon);
Decision_State_startup{NST} = binvar(Number_gen,Number_horizon);
Decision_State_shutdown{NST} = binvar(Number_gen,Number_horizon);
end
%-----------------------------binary variabels----------------------------%
%
%% --------------------------continue variabels---------------------------%
Decision_Generation = cell(Number_scenarios_train,1);
Decision_Generation_cost = cell(Number_scenarios_train,1);
Decision_delta = cell(Number_scenarios_train,1);
Decision_Generation_shedding = cell(Number_scenarios_train,1);
Decision_RES_spilling = cell(Number_scenarios_train,1);
Decision_RES_injected = cell(Number_scenarios_train,1);
Decision_Load_shedding = cell(Number_scenarios_train,1);
Decision_Load_satisfied = cell(Number_scenarios_train,1);
for NST = 1:Number_scenarios_train
Decision_Generation{NST} = sdpvar(Number_gen,Number_horizon,'full');
Decision_Generation_cost{NST} = sdpvar(Number_gen,Number_horizon,'full');
Decision_delta{NST} = sdpvar(Number_gen*(Number_seg+1),Number_horizon,'full');
Decision_Generation_shedding{NST} = sdpvar(Number_gen,Number_horizon,'full');
Decision_RES_spilling{NST} = sdpvar(Number_RES,Number_horizon,'full');
Decision_RES_injected{NST} = sdpvar(Number_RES,Number_horizon,'full');
Decision_Load_shedding{NST} = sdpvar(Number_load,Number_horizon,'full');
Decision_Load_satisfied{NST} = sdpvar(Number_load,Number_horizon,'full');
end
%--------------------------continue variabels-----------------------------%
%
%% ---------------------constraints: general part-------------------------%
Constraint_General = cell(Number_scenarios_train,1);
for NST = 1:Number_scenarios_train
Constraint_General{NST} = [];
end
for NST = 1:Number_scenarios_train
% CC_General 01 - CC_General_09 are general for every model.
% General_01: Non-negative
Constraint_General{NST} = [Constraint_General{NST}, Decision_RES_spilling{NST} >= 0];
Constraint_General{NST} = [Constraint_General{NST}, Decision_RES_injected{NST} >= 0];
Constraint_General{NST} = [Constraint_General{NST}, Decision_Load_shedding{NST} >= 0];
Constraint_General{NST} = [Constraint_General{NST}, Decision_Load_satisfied{NST} >= 0];
%
% General_02: Locally ideal formulation of operation cost
for t = 1:Number_horizon
for i = 1:Number_gen
Constraint_General{NST} = [Constraint_General{NST}, Decision_Generation_cost{NST}(i,t) == Price_split_point(i,:)*Decision_delta{NST}((i-1)*(Number_seg+1)+1:i*(Number_seg+1),t)];
Constraint_General{NST} = [Constraint_General{NST}, sum(Decision_delta{NST}((i-1)*(Number_seg+1)+1:i*(Number_seg+1),t)) == Decision_State_onoff{NST}(i,t)];
Constraint_General{NST} = [Constraint_General{NST}, Generation_split_point(i,1:Number_seg+1)*Decision_delta{NST}((i-1)*(Number_seg+1)+1:i*(Number_seg+1),t) == Decision_Generation{NST}(i,t)];
end
end
Constraint_General{NST} = [Constraint_General{NST}, 0 <= Decision_delta{NST} <= 1];
%
% General_03: Generation limit
for t = 1:Number_horizon
Constraint_General{NST} = [Constraint_General{NST}, Decision_State_onoff{NST}(:,t).*Cus_Data_Gen_capacity(:,3) <= Decision_Generation{NST}(:,t) <= Decision_State_onoff{NST}(:,t).*Cus_Data_Gen_capacity(:,4)];
end
%
% General_04: Startup/shutdown relation with commitment variable
Constraint_General{NST} = [Constraint_General{NST}, Decision_State_startup{NST}(:,1) - Decision_State_shutdown{NST}(:,1) == Decision_State_onoff{NST}(:,1)];
for t = 2:Number_horizon
Constraint_General{NST} = [Constraint_General{NST}, Decision_State_startup{NST}(:,t) - Decision_State_shutdown{NST}(:,t) == Decision_State_onoff{NST}(:,t) - Decision_State_onoff{NST}(:,t-1)];
end
%
% General_05: Min on/off
for i = 1:Number_gen
for t = Cus_Data_Gen_capacity(i,5):Number_horizon
Constraint_General{NST} = [Constraint_General{NST}, sum(Decision_State_startup{NST}(i,t-Cus_Data_Gen_capacity(i,5)+1:t)) <= Decision_State_onoff{NST}(i,t)];
end
end
for i = 1:Number_gen
for t = Cus_Data_Gen_capacity(i,6):Number_horizon
Constraint_General{NST} = [Constraint_General{NST}, sum(Decision_State_shutdown{NST}(i,t-Cus_Data_Gen_capacity(i,6)+1:t)) <= 1 - Decision_State_onoff{NST}(i,t)];
end
end
%
% General_06: Ramping limit
Constraint_General{NST} = [Constraint_General{NST}, Decision_Generation{NST}(:,1) <= Cus_Data_Gen_capacity(:,9).*Decision_State_onoff{NST}(:,1) + Cus_Data_Gen_capacity(:,4).*(1-Decision_State_onoff{NST}(:,1))];
for t = 2:Number_horizon
Constraint_General{NST} = [Constraint_General{NST}, Decision_Generation{NST}(:,t) - Decision_Generation{NST}(:,t-1) <= Cus_Data_Gen_capacity(:,7).*Decision_State_onoff{NST}(:,t-1) + Cus_Data_Gen_capacity(:,9) .*(Decision_State_onoff{NST}(:,t) - Decision_State_onoff{NST}(:,t-1)) + Cus_Data_Gen_capacity(:,4).*(1-Decision_State_onoff{NST}(:,t))];
Constraint_General{NST} = [Constraint_General{NST}, Decision_Generation{NST}(:,t-1)- Decision_Generation{NST}(:,t) <= Cus_Data_Gen_capacity(:,8).*Decision_State_onoff{NST}(:,t) + Cus_Data_Gen_capacity(:,10).*(Decision_State_onoff{NST}(:,t-1) - Decision_State_onoff{NST}(:,t)) + Cus_Data_Gen_capacity(:,4).*(1-Decision_State_onoff{NST}(:,t-1))];
end
%
% General_07: Power balance
for t = 1:Number_horizon
Constraint_General{NST} = [Constraint_General{NST}, sum(Decision_Generation{NST}(:,t) - Decision_Generation_shedding{NST}(:,t)) + sum(Decision_RES_injected{NST}(:,t)) == sum(Decision_Load_satisfied{NST}(:,t))];
end
%
% General_08: Line capacity limit
for t = 1:Number_horizon
Constraint_General{NST} = [Constraint_General{NST}, -Cus_Data_Branch(:,4) <= Cus_Data_Gen_PTDF*(Decision_Generation{NST}(:,t) - Decision_Generation_shedding{NST}(:,t)) + Cus_Data_RES_PTDF*Decision_RES_injected{NST}(:,t) - Cus_Data_Load_PTDF*Decision_Load_satisfied{NST}(:,t) <= Cus_Data_Branch(:,4)];
end
%
% General_09: Generation shedding limits
for t = 1:Number_horizon
Constraint_General{NST} = [Constraint_General{NST}, 0 <= Decision_Generation_shedding{NST}(:,t) <= Decision_Generation{NST}(:,t)];
end
%
end
%% ---------------------constraints: special part-------------------------%
Constraint_Special = cell(Number_scenarios_train,1);
% CC_Special_01 - CC_Special_02 are the special part for specifical model.
for NST = 1:Number_scenarios_train
Constraint_Special{NST} = [];
end
for NST = 1:Number_scenarios_train
% CC_Special_01: RES spilling limits
% CC_Specoal_02: Load shedding limits
% RES_temp = Cus_Data_RES_CUM_00((NST-1)*Number_horizon+1:Number_horizon*NST,5:6)';
% Constraint_Special{NST} = [Constraint_Special{NST}, Decision_RES_spilling{NST}(:) + Decision_RES_injected{NST}(:) == RES_temp(:)];
Constraint_Special{NST} = [Constraint_Special{NST}, Decision_RES_spilling{NST}(:) + Decision_RES_injected{NST}(:) == B(1:Number_horizon*Number_RES,:)*Feature_D_train(:,NST)];
for t = 1:Number_horizon
Constraint_Special{NST} = [Constraint_Special{NST}, Decision_Load_shedding{NST}(:,t) + Decision_Load_satisfied{NST}(:,t) == Realization_load_bus_train(Number_horizon*(NST-1)+t,:)'];
end
end
%% ----------------------------objective----------------------------------%
for NST = 1:Number_scenarios_train
Constraint = Constraint + Constraint_General{NST} + Constraint_Special{NST};
end
Obj_cost_startup = sdpvar(Number_scenarios_train,1);
Obj_cost_shutdown = sdpvar(Number_scenarios_train,1);
Obj_cost_generation = sdpvar(Number_scenarios_train,1);
Obj_cost_generation_shedding = sdpvar(Number_scenarios_train,1);
Obj_cost_RES_spilling = sdpvar(Number_scenarios_train,1);
Obj_cost_load_shedding = sdpvar(Number_scenarios_train,1);
Obj_cost_all = sdpvar(Number_scenarios_train,1);
Obj_loss_abs = sdpvar(Number_scenarios_train,1);
Obj_ERM = lamda*norm(B,1);
for NST = 1:Number_scenarios_train
Constraint = Constraint + [Obj_cost_startup(NST,1) == Cus_Data_Gen_price(:,5).*sum(sum(Decision_State_startup{NST}))];
Constraint = Constraint + [Obj_cost_shutdown(NST,1) == Cus_Data_Gen_price(:,6).*sum(sum(Decision_State_shutdown{NST}))];
Constraint = Constraint + [Obj_cost_generation(NST,1) == sum(Decision_Generation_cost{NST}(:))];
Constraint = Constraint + [Obj_cost_generation_shedding(NST,1) == sum(GS_price*Decision_Generation_shedding{NST}(:))];
Constraint = Constraint + [Obj_cost_RES_spilling(NST,1) == sum(RESS_price*Decision_RES_spilling{NST}(:))];
Constraint = Constraint + [Obj_cost_load_shedding(NST,1) == sum(LS_price*Decision_Load_shedding{NST}(:))];
Constraint = Constraint + [Obj_cost_all(NST,1) == Obj_cost_startup(NST,1) + Obj_cost_shutdown(NST,1) + Obj_cost_generation(NST,1)...
+ Obj_cost_generation_shedding(NST,1) + Obj_cost_RES_spilling(NST,1) + Obj_cost_load_shedding(NST,1)];
Constraint = Constraint + [Obj_loss_abs(NST,1) == Obj_cost_all(NST,1) - z(NST,1)] + [Obj_cost_all(NST,1) >= z(NST,1)];
end
Obj_ERM = Obj_ERM + sum(Obj_loss_abs)/Number_scenarios_train;
%
%% ----------------------------solve it-----------------------------------%
% Warm start
for NST = 1:Number_scenarios_train
assign(Decision_State_onoff{NST}, Recorder_S_onoff{NST});
assign(Decision_State_startup{NST}, Recorder_S_su{NST});
assign(Decision_State_shutdown{NST},Recorder_S_sd{NST});
end
ops = sdpsettings('solver', 'Gurobi','usex0',1);
ops.gurobi.MIPGap = 0.05;
optimize(Constraint,Obj_ERM,ops)
Thanks,
Xianbang
-
Official comment
This post is more than three years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum. Or why not try our AI Gurobot?.
Post is closed for comments.
Comments
1 comment