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 程序结果
