基于MATLAB的绝热压缩空气储能系统代码咨询

Question

想要一份基于MATLAB的绝热压缩空气储能系统容量配置分析的代码，并用NSGA-Ⅱ算法为引擎追求其最优 Pareto 解
储能系统部分代码为（有错误，不会改）

 % Ces_rata--系统额定容量 Pes_rata--系统额定功率
% sp--电网分时电价 sw--上网电价 PW--风机输出功率 PL--用户负载需求功率
% 基本参数
Ces_rata = x(1);
Pes_rata = x(2);
n = 3; % 风机数量
ICpw = 830; % 风电投资成本
PWr = 1500; % 风机额定功率
T = 24; % 小时数
r =0.08;  % 投资利率
per = 1730.97;% CO2排放率 kg/ton
fc = 0.0034;% CO2 排放标准 $/kg 
OMc =200; % CAES运行成本 每年 MWh $
OMw = 0.0122; % 风电运行成本 每年 kwh $
Nday = 230; % 一年内系统运行天数
L = 20; % CAES 的寿命
yta_ch = 0.83;% 压缩部分等熵效率
yta_dis = 0.85;% 膨胀部分等熵效率
yta_1 = 0.1; % CAES 容量下限
yta_2 = 0.9; % CAES 容量上限
Invp = 700;% CAES初始投资成本 KW $
Invc = 5; % CAES初始投资成本 KWh $
SOC_B = [yta_1,yta_2];% SOC储气容量存储状态
SOC(1)=0.5; % SOC初始状态
pa = 0.1; % 大气压强 MPa
ps_max = 6.0; % 储气罐最大压力状态
ps_min = 4.0; % 储气罐最小压力状态
delta_t_soc = 1; % 时间差
Rg = 287.1; % 气体常数
%% 循环遍历
for i = 2:1:T+1 % 历遍小时循环
 delta_P(i-1) = PW(i-1) - PL(i-1); 
 if delta_P(i-1)>0 %先给储能场充电
 soc_a0 = SOC(i-1);%上一时刻 SOC
 soc_1(i-1) = soc_a0 + yta_ch*delta_P(i-1)*delta_t_soc/Ces_rata; % 理论储能状态
 soc_a1 = max([min([soc_1(i-1),SOC_B(2)]),SOC_B(1)]); % 符合装置约束的储能状态
 Pc1 = (soc_a1 - soc_a0)*Ces_rata/yta_ch/delta_t_soc; %最优充电
 Pc1 = min([max([Pc1,0]),Pes_rata]);% 符合功率限制的充电功率
 soc_a1 = max([min([soc_a0 + yta_ch*Pc1*delta_t_soc/Ces_rata,SOC_B(2)]),SOC_B(1)]); % t时刻储气罐储存状态
 Pch(i-1) = Pc1 ;
 Pdis(i-1)=0;
 uch(i-1)=1;% 充电变量
 delta_P1 = delta_P(i-1) - Pc1;
 if delta_P1>0
 Pdu(i-1) = delta_P1 ;% 储能后弃风电量
 Pg(i-1)=0;
end 
elseif delta_P(i-1) < 0 % 储能系统放电
delta_P(i-1) = abs(delta_P(i-1)); % 取绝对值
soc_a0 = SOC(i-1); 
soc_1(i-1) = soc_a0-delta_P(i-1)*delta_t_soc/(yta_dis*Ces_rata); % 理论储能状态
soc_a1=min([max([soc_1(i-1),SOC_B(1)]),SOC_B(2)]); 
Pc1 = (-1*(soc_a1 - soc_a0)*Ces_rata*yta_dis)/delta_t_soc; %最优的放电电量 
Pc1 = min([max([Pc1,0]),Pes_rata]);%找出最优的放电功率 
soc_a1 = min([max([soc_a0 - Pc1*delta_t_soc*yta_dis/Ces_rata,SOC_B(1)]),SOC_B(2)]); %满足 CAES 充放电功率约束的 SOC 
Pdis(i-1) = Pc1 ;% 实际放电功率
udis(i-1) = 1;% 放电变量
Pc1 = abs(Pc1); 
delta_P1 = delta_P(i-1) - Pc1; 
if delta_P1>0 
Pg(i-1) = delta_P1 ; % CAES放电后仍然缺失供给的用户载荷数
end 
end 
SOC(i) = soc_a1; 
Pyh(i-1) = PL(i-1)-Pg(i-1)- Pdis(i-1); % 风电场满足用户侧负荷需求功率
end 
ps = mean(SOC)*(ps_max-ps_min)+ps_min; % 储气罐压力状态
I1 = Nday*sum(Pdis.*sp); % CAES充放电运行效益
I2 = Nday*sum(Pyh.*sw); % 风电上网收益
Qr = 350*sum(Pyh+Pdis)*per/1000;  % 典型日污染物排放量 kg
I3 = Nday*Qr*fc; % 减少排放收益
Icp = 1000*Invp*(r*(1+r)^L)/((1+r)^L-1);  % Pes_rata相关成本系数
Icsoc = 1000*Invc*(r*(1+r)^L)/((1+r)^L-1) ;  % Ces_rata相关成本系数
Acc = Pes_rata*Icp+Ces_rata*Icsoc+ICpw*n*PWr; % 初始投资成本
Aomc = OMc*Ces_rata+n*PWr*OMw; % 运营维护成本
Ag = Nday*sum(Pg.*sp); % 电力不足购电成本
ROI = (I1+I2+I3)/(Acc+Aomc+Ag);  % 投资收益率
Vs = 3600*Ces_rata/(ps*log(ps/pa)); % 储气罐体积

Answer 1

引自免费微信小程序：皆我百晓生

针对你提供的MATLAB代码，首先纠正并优化代码中的一些错误，然后整合NSGA-II算法进行多目标优化。这里，我们不仅修正储能系统部分的代码，还将加入NSGA-II算法的实现，以找到最优的Pareto解。

修正储能系统部分的代码

在你提供的代码中，存在一些逻辑和语法上的小错误。例如，soc_1 和 soc_a1 的计算应更精确地反映系统的物理限制。此外，代码中缺少一些必要的初始化步骤和边界条件检查。下面是修改后的代码段落：

% 初始化变量
delta_P = zeros(1,T);
Pch = zeros(1,T);
Pdis = zeros(1,T);
Pdu = zeros(1,T);
Pg = zeros(1,T);
Pyh = zeros(1,T);
uch = zeros(1,T);
udis = zeros(1,T);
SOC = zeros(1,T+1);

% 循环遍历
for i = 2:1:T+1
    delta_P(i-1) = PW(i-1) - PL(i-1);
    if delta_P(i-1) > 0
        % 充电逻辑
        soc_a0 = SOC(i-1);
        soc_1(i-1) = soc_a0 + yta_ch * delta_P(i-1) * delta_t_soc / Ces_rata;
        soc_a1 = max(min(soc_1(i-1), SOC_B(2)), SOC_B(1));
        Pc1 = (soc_a1 - soc_a0) * Ces_rata / yta_ch / delta_t_soc;
        Pc1 = min(max(Pc1, 0), Pes_rata);
        SOC(i) = max(min(soc_a0 + yta_ch * Pc1 * delta_t_soc / Ces_rata, SOC_B(2)), SOC_B(1));
        Pch(i-1) = Pc1;
        uch(i-1) = 1;
        Pdu(i-1) = max(0, delta_P(i-1) - Pc1);
    elseif delta_P(i-1) < 0
        % 放电逻辑
        delta_P(i-1) = abs(delta_P(i-1));
        soc_a0 = SOC(i-1);
        soc_1(i-1) = soc_a0 - delta_P(i-1) * delta_t_soc / (yta_dis * Ces_rata);
        soc_a1 = min(max(soc_1(i-1), SOC_B(1)), SOC_B(2));
        Pc1 = (-1 * (soc_a1 - soc_a0) * Ces_rata * yta_dis) / delta_t_soc;
        Pc1 = min(max(Pc1, 0), Pes_rata);
        SOC(i) = min(max(soc_a0 - Pc1 * delta_t_soc * yta_dis / Ces_rata, SOC_B(1)), SOC_B(2));
        Pdis(i-1) = abs(Pc1);
        udis(i-1) = 1;
        Pg(i-1) = max(0, delta_P(i-1) - Pc1);
    end
    Pyh(i-1) = PL(i-1) - Pg(i-1) - Pdis(i-1);
end

整合NSGA-II算法

为了使用NSGA-II算法寻找最优的Pareto解，我们将使用MATLAB自带的gamultiobj函数。下面是一个简化的NSGA-II算法实现示例，用于优化储能系统配置：

% NSGA-II 设置
nvars = 2; % 变量数量：Ces_rata 和 Pes_rata
lb = [100, 100]; % 下界
ub = [10000, 10000]; % 上界
options = optimoptions('gamultiobj', 'PopulationSize', 200, 'Generations', 50, 'PlotFcn', @gaplotpareto);

% 目标函数定义
fun = @(x) [computeObjective1(x), computeObjective2(x)]; % 自定义目标函数

% 运行NSGA-II
[solution, ~, ~, ~, ~, ~, ~, ~, ~, ParetoFront] = gamultiobj(fun, nvars, [], [], [], [], lb, ub, [], options);

% 目标函数实现
function F = computeObjective1(x)
    % 这里是第一个目标函数的实现，例如：最小化投资成本
    F = Acc + Aomc + Ag;
end

function F = computeObjective2(x)
    % 这里是第二个目标函数的实现，例如：最大化投资回报率
    F = ROI;
end

注意事项

1. 目标函数：你需要根据实际应用调整computeObjective1和computeObjective2函数，以反映你的具体优化目标。
2. 参数调整：gamultiobj函数中的参数，如PopulationSize和Generations，可以根据需要进行调整以获得更好的优化结果。
3. 边界条件：确保所有变量的边界条件都合理设置，以防止算法陷入不现实的解空间。

通过以上步骤，你应该能够得到一个基于MATLAB的绝热压缩空气储能系统容量配置分析，并使用NSGA-II算法找到最优的Pareto解。