计及调频成本和荷电状态恢复的多储能系统调频功率双层优化【蓄电池经济最优目标下充放电】（基于matlab+yalmip+cplex的蓄电池出力优化）_计及调频成本和荷电状态恢复的 多储能系统调频功率双层优化

摘要：针对电网中不同类型储能电站调频成本、剩余调频
能力存在差异、储能电站内部储能单元 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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146

运行结果如下：