加了一个储能系统后代码跑不了了，求帮助

[dontsleeptonight]

% 清除工作区

clear

clc

%% 参数设置

% 时间向量

T = 24;

% 电力设备参数[光伏PV,风力WT,燃料FC,热电联产CHP]

P_min = [0; 0; 0; 30]; % 最小功率/kW

P_max = [120; 120; 130; 200]; % 最大功率/kW

P_slope_max = [inf; inf; 12; 17; 10]; % 最大爬升率/kWh

efficiency = [NaN; NaN; 0.5; 0.41; 0.9]; % 效率

C_inv = [3000; 6000; 1000; 2000; 600]; % 总投资成本/元

lifetime = [20; 20; 15; 15; 20]; % 寿命/年

lbfc = P_min(3) * ones(T, 1);

ubfc = P_max(3) * ones(T, 1);

lbchp = P_min(4) * ones(T, 1);

ubchp = P_max(4) * ones(T, 1);

% 光伏发电数据2022 年，分布式光伏系统运维成本为0.048 元/ (W·年)=0.048元/(8.76KWh)=0.0055元/度

PV_output = [0; 0; 0; 0; 0; 0; 18.287292818; 56.70349908; 80.39287907; 96.9674647; ...

127.9558011; 146.79558011; 149.00552486; 144.34008594; 122.43093923; ...

100.95150399; 66.08778392; 32.57458564; 14.542664211; 0; 0; 0; 0; 0];%kw

KOMPV=0.0055;

% 风力发电数据我们测算风力发电的度电成本约 0.2 元，考虑期间费用、资产减值以及所得税等因素的完全成本约 0.30-0.35 元/度

WT_output = [153.95948435; 151.56537753; 153.95948435; 157.45856354; 150.64456722; ...

144.5058318; 148.06629834; 136.64825046; 132.0441989; 133.5174954; ...

135.72744015; 131.92142419; 133.08778392; 124.8004911; 129.28176796; ...

121.73112339; 124.1252302; 134.80662983; 144.93554328; 155.61694291; ...

165.74585635; 169.67464702; 170.41129527; 169.24493554];%kw

KOMWT=0.30;

% 购电价和售电价

buy_price = [0.253; 0.234; 0.165; 0.163; 0.2; 0.256; 0.343; 0.406; 0.662; ...

1.02; 2.01; 2.03; 1.58; 1; 0.475; 0.472; 1.02; 1.45; 2.03; ...

2.01; 1.22; 0.824; 0.482; 0.363];

sell_price = [0.18; 0.145; 0.151; 0.151; 0.143; 0.212; 0.252; 0.317; 0.522; ...

0.798; 1.59; 1.59; 1.21; 0.78; 0.401; 0.403; 0.812; 1.13; ...

1.61; 1.59; 0.962; 0.646; 0.392; 0.221];%元/kw

% 环境成本参数[MG,FC,CHP]

CO2_emission_coeff = [890; 441; 670]; % CO2 排放系数，单位 g/(kW·h)

SO2_emission_coeff = [1.8; 0.0022; 0.0036]; % SO2 排放系数，单位 g/(kW·h)

NOx_emission_coeff = [1.6; 0.0136; 0.186]; % NOx 排放系数，单位 g/(kW·h)

CO2_penalty_coeff = 0.21; % CO2 排放惩罚系数，单位 元/kg

SO2_penalty_coeff = 14.84; % SO2 排放惩罚系数，单位 元/kg

NOx_penalty_coeff = 62.96; % NOx 排放惩罚系数，单位 元/kg

% 定义决策变量与其相关参数

P_fc = sdpvar(T, 1); % 燃料电池功率

P_chp = sdpvar(T, 1); % 热电联产组功率

P_buy = sdpvar(T, 1); % 购买电力功率

P_sell = sdpvar(T, 1); % 销售电力功率

buy_indicator = binvar(T, 1); % 购买指示器（二进制变量）

sell_indicator = binvar(T, 1); % 销售指示器（二进制变量）

CCH4 = 4.95; % 天然气价格元/m3

LHVng = 9.72; % 低热值

KOMFC = 0.21; % 燃料电池运维系数

Cng = CCH4 / LHVng; % 天然气价格与其低热值之比

KOMMT = 0.039; % 涡轮机运维系数

%小区电力负荷

P_load = [172.9828042;175.9589947;177.4470899;174.6362434;216.3029101; ...

247.5529101;253.505291;223.5780423;269.8743386;307.0767196;342.6256614;...

311.5410053;260.9457672;210.1851852;216.1375661;259.457672;317.4933862;...

341.3029101;330.7208995;309.8875661;266.8981481;222.255291;196.792328;180.231211];

%% 定义目标函数

% 投资成本

investment_cost = sum(C_inv);

% 风力和光伏运维成本

wt_pv_om_cost = sum(KOMWT .* WT_output) + sum(KOMPV .* PV_output);

% 燃料电池和热电联产成本计算

fc_cost = CCH4/(efficiency(3)*LHVng)*sum(P_fc) + KOMFC*sum(P_fc) ;

chp_cost =Cng/efficiency(4)*sum(P_chp) + KOMMT*sum(P_chp);

% 环境成本

emission_cost =sum(CO2_emission_coeff(1)*CO2_penalty_coeff/1000*(P_buy+P_sell));%环境成本

emission_cost =emission_cost + sum(SO2_emission_coeff(1)*SO2_penalty_coeff/1000*(P_buy+P_sell));

emission_cost =emission_cost + sum(NOx_emission_coeff(1)*NOx_penalty_coeff/1000*(P_buy+P_sell));

emission_cost =emission_cost + sum(CO2_emission_coeff(2)*CO2_penalty_coeff/1000*(P_fc));

emission_cost =emission_cost + sum(SO2_emission_coeff(2)*SO2_penalty_coeff/1000*(P_fc));

emission_cost =emission_cost + sum(NOx_emission_coeff(2)*NOx_penalty_coeff/1000*(P_fc));

emission_cost =emission_cost + sum(CO2_emission_coeff(3)*CO2_penalty_coeff/1000*(P_chp));

emission_cost =emission_cost + sum(SO2_emission_coeff(3)*SO2_penalty_coeff/1000*(P_chp));

emission_cost =emission_cost + sum(NOx_emission_coeff(3)*NOx_penalty_coeff/1000*(P_chp));

%购电成本

purch =sum(P_buy.*buy_price-P_sell.*sell_price);

% 总成本计算，将上述各项成本相加

total_cost = investment_cost + wt_pv_om_cost + fc_cost + chp_cost + emission_cost + purch;

% 定义约束条件

cons=[];

% 1，定义功率平衡约束

for t = 1:T

cons = [cons;P_buy(t) - P_sell(t) + PV_output(t) + WT_output(t) + P_fc(t) + P_chp(t) == P_load(t)];

end

% 2，燃料电池，热电联产组和购售电功率的约束

for t = 1:T

cons = [cons, P_fc(t) >= lbfc(t)];

cons = [cons, P_fc(t) <= ubfc(t)];

cons = [cons, P_chp(t) >= lbchp(t)];

cons = [cons, P_chp(t) <= ubchp(t)];

cons = [cons, P_buy(t) >= 0];

cons = [cons, P_sell(t) >= 0];

end

% 3,燃料电池和热电联产组的功率爬升率约束

if T > 1

% 计算功率差

diff_P_fc = diff(P_fc);

diff_P_chp = diff(P_chp);

% 创建比较用的列向量

slope_limit_fc = P_slope_max(3) * ones(T-1, 1);

slope_limit_chp = P_slope_max(4) * ones(T-1, 1);

% 添加爬升率约束到 cons 数组

cons = [cons; diff_P_fc <= slope_limit_fc];

cons = [cons; -diff_P_fc <= slope_limit_fc]; % 如果需要包括下降斜率

cons = [cons; diff_P_chp <= slope_limit_chp];

cons = [cons; -diff_P_chp <= slope_limit_chp]; % 如果需要包括下降斜率

end

% % 4，定义销售和购买电力的互斥约束

% for t = 1:T

% cons = [cons; buy_indicator(t) + sell_indicator(t) <= 1]; % 购买和销售指示器不能同时为1

% cons = [cons; P_buy(t) <= P_load(t) * buy_indicator(t); % 如果buy_indicator(t)为0，则P_buy(t)必须为0

% P_sell(t) <= sum(P_max) * sell_indicator(t); % 如果sell_indicator(t)为0，则P_sell(t)必须为0

% ];

% end

%%开始求解

%选择求解器

ops = sdpsettings('verbose', 0, 'solver', 'cplex', 'debug', 1);

%调用cplex求解

optimize(cons,total_cost,ops)

% 提取结果

result1 = value(P_fc);

result2 = value(P_chp);

result3 = value(P_buy);

result4 = value(P_sell);

result5 = value(total_cost);

result6 = value(P_load);

result7 = value(P_fc+P_chp+P_buy+WT_output+PV_output);

disp(result5);

% 创建一个新的图形窗口

figure;

% 绘制 P_fc 变量的线状图

plot(result1, 'LineWidth', 1.5); hold on;

% 绘制 P_chp 变量的线状图

plot(result2, 'LineWidth', 1.5);

% 绘制 P_buy 变量的线状图

plot(result3, 'LineWidth', 1.5);

% 绘制 P_sell 变量的线状图

plot(result4, 'LineWidth', 1.5);

% 绘制 P_load 变量的线状图

plot(result6, 'LineWidth', 1.5);

% 绘制 P_load 变量的线状图

plot(result7, 'LineWidth', 1.5);

% 设置图例

legend('P\_fc', 'P\_chp', 'P\_buy', 'P\_sell', 'P\_load','P\_LCG');

% 添加标题和轴标签

title('变量数值线状图');

xlabel('时间（小时）');

ylabel('功率（kW）/ 成本（元）');

% 显示图形

hold off;

 

 

以上是我原本的代码在增加了储能系统部分的建模（如下）后就运行不出答案了

 

Pbat=sdpvar(1,T,'full');%蓄电池出力

Temp_cha=binvar(1,T,'full'); %充电标志

Temp_dis=binvar(1,T,'full'); %放电标志

Temp_static=binvar(1,T,'full'); %电池静置标志

Pcha=sdpvar(1,T);%蓄电池充电电量

Pdis=sdpvar(1,T);%蓄电池放电电量

% 1，定义功率平衡约束（更改后）

for t = 1:T

cons = [cons;P_buy(1,t) - P_sell(1,t) + PV_output(t) + WT_output(t) + P_fc(1,t) + P_chp(1,t) +Pbat(1,t) == P_load(t)];

end

% 4,蓄电池约束

for t=1:24

cons = [cons; P_min(5) <= Pbat(1,t) <= P_max(5)]; % 电池功率范围约束

cons = [cons; 0 <= Pdis(1,t) <= P_max(5)]; % 放电功率非负约束

cons = [cons; P_min(5) <= Pcha(1,t) <= 0]; % 充电功率非正约束

cons = [cons; implies(Temp_cha(1,t),[Pbat(1,t)<=0,Pcha(1,t)==Pbat(1,t),Pdis(1,t)==0])];%充电情况约束

cons = [cons; implies(Temp_dis(1,t),[Pbat(1,t)>=0,Pdis(1,t)==Pbat(1,t),Pcha(1,t)==0])];%放电情况约束

cons = [cons;implies(Temp_static(1,t),[Pbat(1,t)==0,Pdis(1,t)==0,Pcha(1,t)==0])];%静置情况约束

cons = [cons;Temp_cha(1,t)+Temp_dis(1,t)+Temp_static(1,t)==1];

cons = [cons;-170<=sum(Pdis(1,1:t)+Pcha(1,1:t))<=130] ;%SOC约束，电池容量500kwh，初始S0C为54%.20%~80%

end

cons = [cons;sum(Pdis+Pcha)==0] ;%ST=S0，始末SOC相等约束

求大佬帮忙看一下，解释出现问题的原因