目录
1 主要内容
该程序参考《基于全寿命周期成本的配电网蓄电池储能系统的优化配置》全寿命模型和《含分布式发电的微电网中储能装置容量优化配置_刘舒》容量配置部分，主要内容：解决微网中蓄电池优化配置的问题，其中储能的电池容量未知，在一定的框架下对其进行优化，得出满足经济效益最佳的储能容量配置结果，并且在微网的框架下，给出了不同时段的容量配置结果，还给出了微网购电/售电策略，电池充电/放电策略，以及微网中其他单元的配置结果。该程序可以选择采用MATLAB+GUROBI平台求解或者是MATLAB自行求解，对于未安装gurobi的同学更加友好。
程序难点：
该程序主要是采用紧凑约束形式，如Ax<=b/Aeq*x=beq。
这种形式虽然简单，但是理解x具体代表哪些变量还是有些难度。
有个窍门是通过x的输出部分来进一步了解，通过反推方式在理解上述公式的难度就大大降低。x结果输出部分代码如下：
Pbuy_interval=x(1:num_of_hours); Psell_interval=x(num_of_hours+1:2*num_of_hours); Pbattery_inter_disch=x(2*num_of_hours+1:3*num_of_hours); Pbattery_inter_charge=x(3*num_of_hours+1:4*num_of_hours); Px1_interval=x(4*num_of_hours+1:5*num_of_hours); Px2_interval=x(5*num_of_hours+1:6*num_of_hours); Py11_interval=x(6*num_of_hours+1:7*num_of_hours); Py12_interval=x(7*num_of_hours+1:8*num_of_hours); Py13_interval=x(8*num_of_hours+1:9*num_of_hours); Py14_interval=x(9*num_of_hours+1:10*num_of_hours); Py21_interval=x(10*num_of_hours+1:11*num_of_hours); Py22_interval=x(11*num_of_hours+1:12*num_of_hours); Py23_interval=x(12*num_of_hours+1:13*num_of_hours); Py24_interval=x(13*num_of_hours+1:14*num_of_hours); B1_interval=x(14*num_of_hours+1:15*num_of_hours); B2_interval=x(15*num_of_hours+1:16*num_of_hours); B3_interval=x(16*num_of_hours+1:17*num_of_hours); B4_interval=x(17*num_of_hours+1:18*num_of_hours); Z1_interval=x(18*num_of_hours+1:19*num_of_hours); Z2_interval=x(19*num_of_hours+1:20*num_of_hours); Z3_interval=x(20*num_of_hours+1:21*num_of_hours); Z4_interval=x(21*num_of_hours+1:22*num_of_hours);
2 部分程序
read_time=toc %send to display the time it took to read the data from the .xlsx file
num_descrit=length(efficiencies);
%end of input data
%% *****************************************************************************
%prepare the matrixes
icb=initial_capacity_battery;
mincb=minimum_capacity_battery;
maxcb=maximum_capacity_battery;
num_of_hours=length(c);
%1D zero element matrix (time intervals)
zeros_1D=zeros(1,num_of_hours);
%2D zero element matrix (time intervals*time intervals)
zeros_2D=zeros(num_of_hours,num_of_hours);
%Create the positive diagonal matrix with 1.
v = ones(1,num_of_hours);
Diag1_pos = diag(v);
%Create the negative diagonal matrix with -1.
v = ones(1,num_of_hours);
v=v.*-1;
Diag1_neg =diag(v);
%Create the negative diagonal battery efficiency.
Diag_neg_disch_eff=Diag1_neg*battery_effic_disch; %is not implemented yet
Diag_pos_charge_eff=Diag1_pos*(1/battery_effic_charge);
%Create the diagonal matrix OF THE first efficiency
v = ones(1,num_of_hours);
v=v.*efficiencies(1);
Diag_eff1_pos =diag(v);
Diag_eff1_neg=-Diag_eff1_pos;
%Create the diagonal matrix OF THE second efficiency
v = ones(1,num_of_hours);
v=v.*efficiencies(2);
Diag_eff2_pos =diag(v);
Diag_eff2_neg=-Diag_eff2_pos;
%Create the diagonal matrix OF THE third efficiency
v = ones(1,num_of_hours);
v=v.*efficiencies(3);
Diag_eff3_pos =diag(v);
Diag_eff3_neg=-Diag_eff3_pos;
%Create the diagonal matrix OF THE fourth efficiency
v = ones(1,num_of_hours);
v=v.*efficiencies(4);
Diag_eff4_pos =diag(v);
Diag_eff4_neg=-Diag_eff4_pos;
%create lower triangular matrix positive and negative
v=tril(ones(num_of_hours,num_of_hours),-1);
pos_triangular=v+Diag1_pos;
neg_triangular=-pos_triangular;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pos_triangular=pos_triangular/4; %for quarter hour
neg_triangular=neg_triangular/4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Create the diagonal matrix OF THE bounds efficiency
v = ones(1,num_of_hours);
v=v.*bounds_efficiency(1);
Diag_bounds_eff1 =diag(v);
%Create the diagonal matrix OF THE bounds efficiency
v = ones(1,num_of_hours);
v=v.*bounds_efficiency(2);
Diag_bounds_eff2 =diag(v);
%Create the diagonal matrix OF THE bounds efficiency
v = ones(1,num_of_hours);
v=v.*bounds_efficiency(3);
Diag_bounds_eff3 =diag(v);
%Create the diagonal matrix OF THE bounds efficiency
v = ones(1,num_of_hours);
v=v.*bounds_efficiency(4);
Diag_bounds_eff4 =diag(v);
%create the upper bound value of y1 and y2
upper_y1=max_batt_discharge+max(PV);
upper_y2=max_batt_charge*(1/battery_effic_charge);
%Create the diagonal matrix of upper bound value of y1
v = ones(1,num_of_hours);
v=v.*upper_y1; %very big value for ours values
Diag_upper_pos_y1 =diag(v);
Diag_upper_neg_y1 =-Diag_upper_pos_y1;
%Create the diagonal matrix of upper bound value of y2
v = ones(1,num_of_hours);
v=v.*upper_y2; %very big value for ours values
Diag_upper_pos_y2 =diag(v);
Diag_upper_neg_y2 =-Diag_upper_pos_y2;
%% *****************************************************************************
%Objective Function
%
% F=[pgrid_pos pgrid_neg bat PDC y1:y4 B1:B4];
F=[c k*(-1) zeros(1,num_of_hours*4) zeros(1,num_of_hours*4*num_descrit)];
intcon = 6*num_of_hours+num_of_hours*2*num_descrit+1:length(F); %binary variables
% %BOUNDS
%**********Ina***
PgridMax=inf;%limit the power from grid: INA added to avoid spikes in charging batery
%***********
lb=[zeros(1,2*num_of_hours) zeros(1,2*num_of_hours) zeros(1,2*num_of_hours) zeros(1,num_of_hours*4*num_descrit) ];
% ub=[ones(1,2*num_of_hours)*(inf) ones(1,num_of_hours)*(max_batt_discharge) ones(1,num_of_hours)*(max_batt_charge) ones(1,2*num_of_hours)*(inf) ones(1,num_of_hours*4*num_descrit)*(inf) ];
ub=[ones(1,2*num_of_hours)*(PgridMax) ones(1,num_of_hours)*(max_batt_discharge) ones(1,num_of_hours)*(max_batt_charge) ones(1,2*num_of_hours)*(PgridMax) ones(1,num_of_hours*4*num_descrit)*(PgridMax) ];
%Equalities Constraints
Aeq=[Diag1_pos,Diag1_neg,zeros_2D,zeros_2D,zeros_2D,Diag1_neg,Diag1_pos,Diag1_pos,Diag1_pos,Diag1_pos,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D;
zeros_2D,zeros_2D,Diag_neg_disch_eff,Diag_pos_charge_eff,Diag1_pos,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,Diag1_neg,Diag1_neg,Diag1_neg,Diag1_neg,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D ;
zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,zeros_2D,Diag1_pos,Diag1_pos,Diag1_pos,Diag1_pos,Diag1_pos,Diag1_pos,Diag1_pos,Diag1_pos;
];
beq=[Load,PV,ones(1,num_of_hours)];
%Inequalities Constraints
3 程序结果
1008
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
