摘要：针对电网中不同类型储能电站调频成本、剩余调频
能力存在差异、储能电站内部储能单元 SOC 过高或过低
的问题，提出计及调频成本和 SOC 恢复的多储能系统调
频功率双层优化策略，该策略包含调频功率优化层和
SOC 优化层：在调频功率优化层引入与储能电站剩余调
频能力相关的抗拒系数，与调频成本共同构成目标函数，
根据不同类型储能电站的调频成本和剩余调频能力实现
调频功率优化分配；在 SOC 优化层引入 SOC 权重以恢复
储能单元 SOC 为目标对调频功率优化结果进行再分配。
采用包含三种不同电池类型的区域电网验证所提分配策
略的有效性，结果表明：与 3 种对比策略相比，所提策略
不仅具有较好的经济性，还能恢复各储能单元的 SOC，
显著提高了调频效果。

目标函数经济成本最小如下：

目标函数参数：

蓄电池参数：

蓄电池soc初始值：

四种策略方案：

对策略二基于matlab+yalmip+cplex的蓄电池出力优化复现：
蓄电池经济最优目标下充放电
程序如下：

// matlab+yalmip+cplex的蓄电池出力优化
%% 策略二，按经济最优充放电
clc
clear
close all
[agc]=xlsread('agcdata','Sheet1','c2:c181');%agc参数
canshu;
X = sdpvar(8,180);         %电站1
Y = sdpvar(8,180);         %电站2
Z = sdpvar(8,180);         %电站3

%% 约束矩阵
constraint=[];
% objective = 0;

%% 功率平衡约束

    for j=1:180
        if agc(j)>=0
        constraint=[constraint,agc(j) == sum(X(:,j))*e1 + sum(Y(:,j))*e2 + sum(Z(:,j))*e3 ];
        else
        constraint=[constraint,agc(j) == sum(X(:,j))/e1 + sum(Y(:,j))/e2 + sum(Z(:,j))/e3 ];
        end
    end

 %% 储能约束
 for i=1:180
     % 第一单元
     constraint=[constraint,S1*soc_min1 <= S1*soc(1,1) - sum(X(1,1:i)) <= S1*soc_max1];
     constraint=[constraint,S1*soc_min1 <= S1*soc(1,2) - sum(X(2,1:i)) <= S1*soc_max1];
     constraint=[constraint,S1*soc_min1 <= S1*soc(1,3) - sum(X(3,1:i)) <= S1*soc_max1];
     constraint=[constraint,S1*soc_min1 <= S1*soc(1,4) - sum(X(4,1:i)) <= S1*soc_max1];
     constraint=[constraint,S1*soc_min1 <= S1*soc(1,5) - sum(X(5,1:i)) <= S1*soc_max1];
     constraint=[constraint,S1*soc_min1 <= S1*soc(1,6) - sum(X(6,1:i)) <= S1*soc_max1];
     constraint=[constraint,S1*soc_min1 <= S1*soc(1,7) - sum(X(7,1:i)) <= S1*soc_max1];
     constraint=[constraint,S1*soc_min1 <= S1*soc(1,8) - sum(X(8,1:i)) <= S1*soc_max1];
     
     % 第二单元
     constraint=[constraint,S2*soc_min2 <= S2*soc(2,1) - sum(Y(1,1:i)) <= S2*soc_max2];
      constraint=[constraint,S2*soc_min2 <= S2*soc(2,2) - sum(Y(2,1:i)) <= S2*soc_max2];
      constraint=[constraint,S2*soc_min2 <= S2*soc(2,3) - sum(Y(3,1:i)) <= S2*soc_max2];
     constraint=[constraint,S2*soc_min2 <= S2*soc(2,4) - sum(Y(4,1:i)) <= S2*soc_max2]; 
      constraint=[constraint,S2*soc_min2 <= S2*soc(2,5) - sum(Y(5,1:i)) <= S2*soc_max2];
      constraint=[constraint,S2*soc_min2 <= S2*soc(2,6) - sum(Y(6,1:i)) <= S2*soc_max2];
      constraint=[constraint,S2*soc_min2 <= S2*soc(2,7) - sum(Y(7,1:i)) <= S2*soc_max2];
      constraint=[constraint,S2*soc_min2 <= S2*soc(2,8) - sum(Y(8,1:i)) <= S2*soc_max2];
      
      % 第三单元
     constraint=[constraint,S3*soc_min3 <= S3*soc(3,1) - sum(Z(1,1:i)) <= S3*soc_max3];
     constraint=[constraint,S3*soc_min3 <= S3*soc(3,2) - sum(Z(2,1:i)) <= S3*soc_max3];
     constraint=[constraint,S3*soc_min3 <= S3*soc(3,3) - sum(Z(3,1:i)) <= S3*soc_max3];
     constraint=[constraint,S3*soc_min3 <= S3*soc(3,4) - sum(Z(4,1:i)) <= S3*soc_max3];
     constraint=[constraint,S3*soc_min3 <= S3*soc(3,5) - sum(Z(5,1:i)) <= S3*soc_max3];
     constraint=[constraint,S3*soc_min3 <= S3*soc(3,6) - sum(Z(6,1:i)) <= S3*soc_max3];
     constraint=[constraint,S3*soc_min3 <= S3*soc(3,7) - sum(Z(7,1:i)) <= S3*soc_max3];
     constraint=[constraint,S3*soc_min3 <= S3*soc(3,8) - sum(Z(8,1:i)) <= S3*soc_max3];
     
 end
  



%% 成本计算
% 投资成本
inv_cost1=(cost_s1*S1*0.08*1.08.^long1)/(365*1440*((1.08.^long1)-1))  ;
inv_cost2= (cost_s2*S2*0.08*1.08.^long2)/(365*1440*((1.08.^long2)-1)) ;
inv_cost3= (cost_s3*S3*0.08*1.08.^long3)/(365*1440*((1.08.^long3)-1))   ;
obj1=inv_cost1+inv_cost2+inv_cost3;

% 损耗成本
loss_cost=520;
obj2_1=0;
obj2_2=0;
obj2_3=0;
% for i=1:180
%      if agc(i)>=0         
%        obj2_1(i)=loss_cost*(1/e1 - 1)*sum(X(1:8,i))  ;%*REF_1(i)
%        obj2_2(i)=loss_cost*(1/e2 - 1)*sum(Y(1:8,i))  ;
%        obj2_3(i)=loss_cost*(1/e3 - 1)*sum(Z(1:8,i))  ;
%      else
%        obj2_1(i)=loss_cost*(1-e1)*sum(X(1:8,i))  ;%*REF_1(i)
%        obj2_2(i)=loss_cost*(1-e2)*sum(Y(1:8,i))  ;
%        obj2_3(i)=loss_cost*(1-e3)*sum(Z(1:8,i))  ;
%      end
% end

for i=1:180
    for n=1:8
       if agc(i)>=0         
       obj2_1=obj2_1+loss_cost*(1/e1 - 1)*X(n,i)  ;%*REF_1(i)
       obj2_2=obj2_2+loss_cost*(1/e2 - 1)*Y(n,i)  ;
       obj2_3=obj2_3+loss_cost*(1/e3 - 1)*Z(n,i)  ;
       else
       obj2_1=obj2_1+loss_cost*(1-e1)*X(n,i)  ;%*REF_1(i)
       obj2_2=obj2_2+loss_cost*(1-e2)*Y(n,i)  ;
       obj2_3=obj2_3+loss_cost*(1-e3)*Z(n,i)  ;
       end
     end
end
obj2=obj2_1 + obj2_2 + obj2_3;

% 折损成本

life_cost1 = 0.5*cost_p1*S1/times1 ;
life_cost2 = 0.5*cost_p2*S2/times2;
life_cost3 = 0.5*cost_p3*S3/times3 ;
obj3_1=0;
obj3_2=0;
obj3_3=0;
for i=1:180
    for n=1:8
     if agc(i)>=0
        obj3_1=obj3_1+life_cost1*(X(n,i))/(S1*e1*3600) ;
        obj3_2=obj3_2+life_cost2*(Y(n,i))/(S2*e2*3600) ;
        obj3_3=obj3_3+life_cost3*(Z(n,i))/(S1*e3*3600) ;
       else
        obj3_1=obj3_1+life_cost1*(X(n,i))*e1/(S1*3600);
        obj3_2=obj3_2+life_cost1*(X(n,i))*e2/(S1*3600);
        obj3_3=obj3_3+life_cost1*(X(n,i))*e3/(S1*3600);
     end
    end
end
obj3=obj3_1+obj3_2+obj3_3;
% obj4=obj3_1*(REF_1);
%% 目标函数
      objective=obj1+obj2 +obj3;
% objective=sum(objective_);

 %% 求解
ops=sdpsettings('verbose', 1, 'solver', 'cplex');
sol=optimize(constraint,objective,ops);
objective=value(objective)
%% 输出
figure
plot(agc)
X=value(X);
Y= value(Y);
Z=value(Z);
P=zeros(1,180);
for i=1:180
    for n=1:8
        P(i)=P(i)+X(n,i)+Y(n,i)+Z(n,i);
    end
end
hold on
plot(P)

运行结果如下：