I P v 6的分段

IPv6的分段只能由源节点和目的节点进行,这样就简化了包头并减少了用于选路的开销。逐跳分段被认为是一种有害的方法。首先,它在端到端的分段中将产生更多的分段。此外在传输中,一个分段的丢失将导致所有分段重传。IPv6的确可以通过其扩展头来支持分段,但是如下所述,了解IPv4分段如何工作将有助于了解IPv6中为什么要进行改变。

IPv4中,当一个没有分段的包由于太长而无法沿着发送源到目的地的网络链路进行传输时,就需要进行包的分段。举例来说,一个源节点可以创建一个长度为1500字节的包,并把它向Internet上的某个远端目的地发送。这个包通过源节点的本地以太网到达该节点的默认路由器。然后路由器通过其链路把数据发到Internet上,这条链路可能是到一个ISP的点到点连接。在Internet中的某处或离目的节点较近的某处,可能有条网络链路无法处理这样一大块的数据。在这种情况下,使用该网络链路的路由器将不得不把1500字节的数据报分割成许多不超过下一个网络的最大传输单元(MTU)的分段。因此,如果假设下一个链路可以处理的包长度不能超过1280字节的话,路由器将把最初的一个包分割为两个。第一个包的长度为1260字节,留下的20字节用于IPv4头。第二段的长度就是剩余数据的长度,240字节,另外再用20字节作为另一个IPv4头。

IPv4中的分段由包沿途的中间路由器根据需要进行。进行分段的路由器根据需要修改包头并在其中包含进最初的包的数据报标识,同时还将正确地设置分段标志和分段偏移值。当目的节点收到由此产生的分段包之后,该系统必须根据每个分段包的IPv4头后的分段数据重组最初的包。

在使用了分段之后,不论中间的网络是什么类型,不同类型网络上的节点都可以互操作,源节点无需了解任何有关目的节点网络的信息,同时也无需了解它们之间的网络信息。这一直被认为是一个不错的特性,由于不需要节点或路由器存储信息或记录整个Internet的结构,从而Internet可以获得很好的扩展性。但另一方面,它也为路由器带来了性能方面的问题,对IP包进行分段消耗了沿途路由器和目的地的处理能力和时间。了解IP数据报标识、计算分段偏移值、真正把数据分段以及在目的地进行重装都会带来额外的开销。

问题在于对于任何一个指定的路由器,虽然源节点能够了解链路的MTU是多大,但却没有办法事先知道整个路径的MTU。路径MTU是源节点和目的节点之间在不分段时可以沿着该路由穿越任何网络的最大包长。

然而,目前有两种方法可以减少或消除对于分段的需求。第一种方法可用在IPv4中,它使用一种叫做“路径MTU发现”的方法。通过这种方法,路由器可以向目的地发送一个包来报告该路由器上链路的MTU值。如果包到达了一条必须对其进行分段的链路,负责分段的路由器将使用ICMP回送一个报文来指出分段路由器上链路的MTU值。这种过程可以重复进行直到路由器确定路径MTU为止。

另一种减少分段需求的方法是要求所有支持IP的链路必须能够处理一些合理的最小长度的包。换句话说,如果一个链路的MTU超过20字节,那么所有的节点都必须准备产生可观数量的分段包。另一方面,如果能够提出所有网络链路都可以适应的某个合理的长度,并把它设置为允许包长度的绝对最小值,那么就可以消灭分段。

IPv6中实际上同时使用了上面两种方法。在最初的RFC中,IPv6规定每个链路支持的MTU最小为576字节。那么这些包的净荷长度将是536字节,另外40字节用于IPv6头。由于RFC1883发表于1995年,后来产生了很多关于更大的MTU的争论。在Huitema提出的报告(参见《IPv6:新的IP》第2版,Prentice-Hall)中,建议值为1997,SteveDeering则正在促使将MTU值改为1500字节。在最新的于1997年11月发表的Internet草案中,MTU值被设为1280字节。很明显,关注的焦点在于:倡导较短MTU的人希望那些不能支持较长MTU的网络不会被完全丢弃,而倡导较长MTU的人不希望为照顾小部分接近于废弃的网络而使得整个Internet的性能下降。

为了对较短的MTU进行一些弥补,IPv6标准中强烈推荐所有IPv6节点都支持路径MTU发现。路径MTU发现最早出现在RFC1191中,其中使用了分段标志中的“不能分段”来要求中间路由器在发现包太长时返回一个ICMP出错报文。

路径MTU发现的IPv6版本在RFC1981(IPv6的路径MTU发现)中描述。这是对原有的RFC1191的升级,但其中加入了一些改变使之可以工作在IPv6中。其中最重要的是,由于IPv6头中不支持分段,因此也就没有“不能分段”位。正在执行路径MTU发现的节点只是简单地在自己的网络链路上向目的地发送允许的最长包。如果一条中间链路无法处理该长度的包,尝试转发路径MTU发现包的路由器将向源节点回送一个ICMPv6出错报文。然后源节点将发送另一个较小的包。这个过程将一直重复,直到不再收到ICMPv6出错报文为止,然后源节点就可以使用最新的MTU作为路径MTU。

这里需要注意,有一些实例并没有实现路径MTU发现。例如,使用最小IPv6实现来进行远程网络启动的终端只是简单地使用576字节的路径MTU。从源节点到目的节点的IPv6分段,作为一个扩展头来实现。

摘自: http://www.cnpaf.net/Class/IPV6/0532918532776699.html 
检查以下代码,指出错误的地方:%% 日前调度模型 clc clear %% 数据导入 % 调度周期 T = 24; % 风电预测功率 P_wt_pre = [419.043583980000,374.323177900000,392.119523040000,398.681756480000,460.739421520000,405.595056020000,235.598101760000,222.689675320000,260.448460160000,277.386108580000,183.566952380000,120.085357742000,81.9212982000000,127.760610824000,163.730879200000,195.313565580000,209.128520300000,267.555362200000,374.097494120000,407.577665040000,421.392893460000,468.046794740000,494.700020360000,479.292094960000]; % 电负荷 P_Load_e = 0.74*[1209.12680000000,1391.49080000000,1417.72400000000,1431.18080000000,1398.85200000000,1322.37000000000,1410.94800000000,1492.84800000000,1554.33600000000,1390.15800000000,1770.80400000000,1786.93200000000,1512.25200000000,1305.10800000000,1210.39800000000,1228.51400000000,1440.85200000000,1770.42600000000,1969.12800000000,1685.12400000000,1488.74600000000,1445.37400000000,1286.13800000000,1377.11000000000]; % 热负荷 P_Load_h = 0.72*[358.676650700000,358.415415000000,370.788488900000,369.577305000000,344.641167600000,321.889910800000,352.311998500000,326.544656500000,315.744330800000,265.474819400000,232.010437700000,227.021414800000,296.575979700000,292.170468200000,267.079475800000,273.994086500000,288.173414800000,270.458472200000,227.607740900000,345.550477800000,322.246141400000,342.399073100000,347.968033700000,370.337263500000]; % 分时电价 C_buy_e = 100*[0.29 0.29 0.29 0.29 0.29 0.69 0.69 0.69 0.84 0.84 0.84 0.69 0.69 0.69 0.69 0.69 0.69 0.69 0.64 0.64 0.55 0.55 0.29 0.29]; %% 设置参数 P_CHP_max = [150,80]; % CHP机组处理(2台) P_MT_max = [400,455,200]; % 火电机组最大出力 P_MT_min = [200,120,100]; % 火电机组最小出力 R_u = 4*[50,50,25]; % 火电机组上爬坡 a = [0.00048,0.00031,0.0002]; % 成本系数a b = [16.2,17.3,16.6]; % 成本系数b c = [1000,970,700]; % 成本系数c C_eg = [0.91,0.95,0.98]; % 碳排放系数 E_bata = 0.81; % 碳捕集效率 P_yita = [1.05,1.05,1.05]; % 再生塔和压缩机最大工作状态系数 P_lamda = [0.269,0.269,0.269]; % 单位捕碳量能耗 P_D1 = 0.5*[20,20,15]; % 碳捕集固定能耗 P_D_gd = ones(3,24); for i=1:3 P_D_gd(i,:)=P_D1(i).*ones(1,24); end M_MEA = 61.08; % M_MEA的摩尔质量 M_co2 = 44; % 二氧化碳的摩尔质量 C_sita = 0.4; % 再生塔解析量 C_CR = 30; % 醇胺溶液浓度(%) C_rou_R = 1.01; % 醇胺溶液密度 V_CR = [2000,2000,2000]; % 溶液储液装置容量 K_R = 1.17; % 乙醇胺溶剂成本系数 C_fai = 1.5; % 溶剂运行损耗系数 % h区间分段 h_min = [200,120,100]; h_max = [401,456,201]; h_gn = 5; hl_1 = (h_max-h_min)./h_gn; hl_2 = zeros(3,h_gn+1); % k区间分段 k_min = [-400,-455,-200]; k_max = [-199,-119,-99]; k_gn = 5; kl_1 = (k_max-k_min)./k_gn; kl_2 = zeros(3,k_gn+1); for i=1:3 hl_2(i,:)=h_min(i):hl_1(i):h_max(i); end for i=1:3 kl_2(i,:)=k_min(i):kl_1(i):k_max(i); end % 热电比系数 heat_ratio = 0.6 / 0.4; % 燃气轮机热电比 %% 决策变量 % 电力源出力 P_CHP_e = sdpvar(2,24,'full'); % 燃气轮机电出力 P_wt_e = sdpvar(1,24,'full'); % 风电机组出力 P_MT_e = sdpvar(3,24,'full'); % 火电机组出力 P_e_EB = sdpvar(2,24,'full'); % 电锅炉耗电 % 热力源出力 P_CHP_h = sdpvar(2,24,'full'); % 燃气轮机热出力 P_EB_h = sdpvar(2,24,'full'); % 电锅炉热出力 % 天然气输入 P_g_CHP = sdpvar(2,24,'full'); % 天然气需求 % 碳捕集相关 E_G_all = sdpvar(3,24,'full'); % 碳捕集机组产生的总碳排放 E_total_co2 = sdpvar(3,24,'full'); % 机组捕获的总碳排放 E_CG_co2 = sdpvar(3,24,'full'); % 储液装置提供的待捕集二氧化碳量 P_B = sdpvar(3,24,'full'); % 机组运行能耗 P_J = sdpvar(3,24,'full'); % 机组的净出力 V_CA = sdpvar(3,24,'full'); V_FY = sdpvar(3,24,'full'); % 富液体积 V_PY = sdpvar(3,24,'full'); % 贫液体积 P_tra_e = sdpvar(1,24,'full'); % 可转移电负荷 P_cut_e = sdpvar(1,24,'full'); % 可削减电负荷 P_DE_Load_e = sdpvar(1,24,'full'); % 经过过需求响应后的电负荷 H_tra_h = sdpvar(1,24,'full'); % 可转移热负荷 H_cut_h = sdpvar(1,24,'full'); % 可削减热负荷 P_DH_Load_h = sdpvar(1,24,'full'); % 经过过需求响应后的热负荷 % 旋转备用变量 spi_reserve = sdpvar(1,24,'full'); % 正出力分段专用 (三维数组存储) y1 = sdpvar(3, 24, 'full'); h_w = sdpvar(h_gn+1, 24, 3, 'full'); % 权重变量 [分段点×时间×机组] h_z = binvar(h_gn, 24, 3, 'full'); % SOS2二进制变量 [分段×时间×机组] % 负出力分段专用 (三维数组存储) y2 = sdpvar(3, 24, 'full'); k_w = sdpvar(k_gn+1, 24, 3, 'full'); % 权重变量 [分段点×时间×机组] k_z = binvar(k_gn, 24, 3, 'full'); % SOS2二进制变量 [分段×时间×机组] %% 约束条件 Cons = []; % 初始化约束集合 % 1. 机组出力约束 for t = 1:T % 电力平衡 Cons = [Cons, sum(P_CHP_e(:, t)) + sum(P_MT_e(:, t)) + P_wt_e(t) == P_DE_Load_e(t) + sum(P_e_EB(:, t))]; % 热力平衡 Cons = [Cons, sum(P_CHP_h(:, t)) + sum(P_EB_h(:, t)) == P_DH_Load_h(t)]; for i = 1:3 % 三个火电机组 % 碳捕集约束 Cons = [Cons, E_G_all(i, t) == C_eg(i) * P_MT_e(i, t), % 总碳排放 0 <= E_CG_co2(i, t), % 待捕集CO2非负 0 <= P_B(i, t), % 碳捕集能耗非负 E_total_co2(i, t) == E_CG_co2(i, t) + 0.25 * E_bata * C_eg(i) * (y1(i, t) - y2(i, t)), % 捕获总碳量 0 <= E_total_co2(i, t) <= P_yita(i) * E_bata * C_eg(i) * P_MT_max(i), % 捕获量范围 P_B(i, t) == P_lamda(i) * E_total_co2(i, t), % 捕集能耗关系 P_MT_e(i, t) == P_D_gd(i, t) + P_J(i, t) + P_B(i, t), % 碳捕集电厂功率平衡 P_MT_min(i) <= P_MT_e(i, t) <= P_MT_max(i), % 火电机组出力范围 P_MT_min(i) - P_lamda(i)*P_yita(i)*E_bata*C_eg(i)*P_MT_max(i) - P_D_gd(i,t) <= P_J(i, t) <= P_MT_max(i) - P_D_gd(i, t) % 净出力范围 ]; end % 风电约束 Cons = [Cons, 0 <= P_wt_e(t) <= P_wt_pre(t), % 风电出力范围 sum(P_e_EB(:, t)) + P_wt_e(t) <= P_wt_pre(t) % 弃风约束 ]; end % 2. 旋转备用约束 for t = 1:T Cons = [Cons, spi_reserve(t) <= sum(R_u), % 爬坡能力限制 spi_reserve(t) <= sum(P_MT_max) - sum(P_MT_e(:, t)), % 可用容量限制 spi_reserve(t) >= 0.01 * P_DE_Load_e(t), % 最小备用要求 spi_reserve(t) >= 0 % 物理意义要求 ]; end % 3. 分段线性化约束 % 正出力分段 (y1) for i = 1:3 for t = 1:T % 二次函数近似 Cons = [Cons, y1(i, t) == (hl_2(i, :).^2) * h_w(:, t, i)]; % SOS2约束 Cons = [Cons, sum(h_w(:, t, i)) == 1, % 权重和为1 h_w(:, t, i) >= 0, % 权重非负 h_w(1, t, i) <= h_z(1, t, i), % 首端点约束 h_w(end, t, i) <= h_z(end, t, i) % 末端点约束 ]; % 中间点约束 for s = 2:h_gn Cons = [Cons, h_w(s, t, i) <= h_z(s-1, t, i) + h_z(s, t, i) ]; end % 二进制变量约束 Cons = [Cons, sum(h_z(:, t, i)) == 1, % 仅一个区间激活 h_z(:, t, i) >= 0 % 二进制变量非负 ]; end end % 负出力分段 (y2) for i = 1:3 for t = 1:T % 二次函数近似 Cons = [Cons, y2(i, t) == (kl_2(i, :).^2) * k_w(:, t, i)]; % SOS2约束 Cons = [Cons, sum(k_w(:, t, i)) == 1, % 权重和为1 k_w(:, t, i) >= 0, % 权重非负 k_w(1, t, i) <= k_z(1, t, i), % 首端点约束 k_w(end, t, i) <= k_z(end, t, i) % 末端点约束 ]; % 中间点约束 for s = 2:k_gn Cons = [Cons, k_w(s, t, i) <= k_z(s-1, t, i) + k_z(s, t, i) ]; end % 二进制变量约束 Cons = [Cons, sum(k_z(:, t, i)) == 1, % 仅一个区间激活 k_z(:, t, i) >= 0 % 二进制变量非负 ]; end end % 4. 爬坡约束 for t = 1:T-1 for i = 1:3 % 火电机组 Cons = [Cons, -R_u(i) <= P_MT_e(i, t+1) - P_MT_e(i, t) <= R_u(i)]; end for j = 1:2 % 燃气轮机和电锅炉 Cons = [Cons, P_CHP_e(j, t+1) - P_CHP_e(j, t) <= 50, P_e_EB(j, t+1) - P_e_EB(j, t) <= 30 ]; end end % 5. 燃气系统约束 for t = 1:T for j = 1:2 Cons = [Cons, 0 <= P_g_CHP(j, t) <= 130, % 供气限制 0 <= P_CHP_e(j, t) <= P_CHP_max(j), % 发电限制 P_CHP_e(j, t) == 0.98 * P_g_CHP(j, t), % 气-电转换 P_CHP_h(j, t) == heat_ratio * P_CHP_e(j, t), % 热电关系 0 <= P_e_EB(j, t) <= P_wt_pre(t) - P_wt_e(t), % 电锅炉限制 P_EB_h(j, t) == 0.98 * P_e_EB(j, t) % 电热转换 ]; end end % 6. 碳捕集动态约束 V_PY0 = [1000, 1000, 1000]; % 贫液初始值 V_FY0 = [1000, 1000, 1000]; % 富液初始值 for i = 1:3 % 初始条件 (t=1) Cons = [Cons, V_CA(i, 1) == (M_MEA * E_CG_co2(i, 1)) / (M_co2 * C_sita * C_CR * C_rou_R), V_PY(i, 1) == V_PY0(i) + V_CA(i, 1), V_FY(i, 1) == V_FY0(i) - V_CA(i, 1) ]; % 后续时间点 (t=2 to T) for t = 2:T Cons = [Cons, V_CA(i, t) == (M_MEA * E_CG_co2(i, t)) / (M_co2 * C_sita * C_CR * C_rou_R), V_FY(i, t) == V_FY(i, t-1) - V_CA(i, t), V_PY(i, t) == V_PY(i, t-1) + V_CA(i, t) ]; end % 容量约束 (t=1 to T) for t = 1:T Cons = [Cons, 0 <= V_FY(i, t) <= V_CR(i), % 富液容量 0 <= V_PY(i, t) <= V_CR(i) % 贫液容量 ]; end % 周期平衡 Cons = [Cons, V_PY0(i) == V_PY(i, T), V_FY0(i) == V_FY(i, T) ]; end % 7. 需求响应约束 % 电负荷响应 for t = 1:T if (t >= 11 && t <= 12) || (t >= 18 && t <= 20) Cons = [Cons, -0.02*P_Load_e(t) <= P_cut_e(t) <= 0]; else Cons = [Cons, P_cut_e(t) == 0]; end end % 热负荷响应 for t = 1:T if (t >= 1 && t <= 3) || (t >= 11 && t <= 12) Cons = [Cons, -0.02*P_Load_h(t) <= H_cut_h(t) <= 0]; else Cons = [Cons, H_cut_h(t) == 0]; end end % 响应平衡约束 Cons = [Cons, -0.08*P_Load_e <= P_tra_e <= 0.08*P_Load_e, % 电转移范围 sum(P_tra_e) == 0, % 电转移守恒 P_DE_Load_e == P_Load_e + P_tra_e + P_cut_e, % 响应后电负荷 -0.08*P_Load_h <= H_tra_h <= 0.08*P_Load_h, % 热转移范围 sum(H_tra_h) == 0, % 热转移守恒 P_DH_Load_h == P_Load_h + H_tra_h + H_cut_h % 响应后热负荷 ]; %% 目标函数 F = 0; % 火电机组燃料费用 (使用二次成本函数) F_gen = 0; for i = 1:3 for t = 1:T F_gen = F_gen + a(i) * P_MT_e(i,t)^2 + b(i) * P_MT_e(i,t) + c(i); end end % 燃气费用 (假设天然气价格为50元/单位) F_gas = 50 * sum(sum(P_g_CHP)); % 碳捕集相关成本 F_carbon = 0; for i = 1:3 for t = 1:T % 溶剂成本 (基于溶液体积) F_carbon = F_carbon + K_R * C_fai * V_CA(i,t); % 碳交易成本 (假设碳交易价格为150元/吨) F_carbon = F_carbon + 150 * (E_G_all(i,t) - E_total_co2(i,t)); end end % 需求响应成本 F_DR_e = 0; % 电能需求响应成本 F_DR_h = 0; % 热能需求响应成本 for t = 1:T F_DR_e = F_DR_e + 0.25 * C_buy_e(t) * abs(P_tra_e(t)) + 0.65 * C_buy_e(t) * abs(P_cut_e(t)); F_DR_h = F_DR_h + 0.25 * abs(H_tra_h(t)) + 0.65 * abs(H_cut_h(t)); end F_DR_all = F_DR_e + F_DR_h; % 总目标函数 F = F_gen + F_gas + F_carbon + F_DR_all; %% CPLEX求解器配置 ops = sdpsettings('solver','cplex','verbose',2,'usex0',0); ops.cplex.mip.tolerances.mipgap = 1e-6; result = optimize(Cons,F,ops); F = value(F); display(['最优规划结果 : ', num2str(F)]);
11-15
检查和优化重构以下代码:%% 日前调度优化模型 clc; clear; %% 参数初始化 T = 24; % 调度周期 % 预测数据 P_wt_pre = 1.4 * [299.3168, 267.3737, 280.0854, 284.7727, 329.0996, 289.7108, 168.2844, 159.0641, ... 186.0346, 198.1329, 131.1193, 85.7753, 58.5152, 91.2576, 116.9506, 139.5097, ... 149.3775, 191.1110, 267.2125, 291.1269, 300.9949, 334.3191, 353.3572, 342.3515]; P_Load_e = 0.74 * [1209.1268, 1391.4908, 1417.7240, 1431.1808, 1398.8520, 1322.3700, 1410.9480, 1492.8480, ... 1554.3360, 1390.1580, 1770.8040, 1786.9320, 1512.2520, 1305.1080, 1210.3980, 1228.5140, ... 1440.8520, 1770.4260, 1969.1280, 1685.1240, 1488.7460, 1445.3740, 1286.1380, 1377.1100]; P_Load_h = 0.7 * [358.6767, 358.4154, 370.7885, 369.5773, 344.6412, 321.8899, 352.3120, 326.5447, ... 315.7443, 265.4748, 232.0104, 227.0214, 296.5760, 292.1705, 267.0795, 273.9941, ... 288.1734, 270.4585, 227.6077, 345.5505, 322.2461, 342.3991, 347.9680, 370.3373]; gama_elec = 100 * [0.29, 0.29, 0.29, 0.29, 0.29, 0.69, 0.69, 0.69, 0.84, 0.84, 0.84, 0.69, ... 0.69, 0.69, 0.69, 0.69, 0.69, 0.69, 0.64, 0.64, 0.55, 0.55, 0.29, 0.29]; % 机组参数 a = [0.00048, 0.00031, 0.0002]; % 成本系数a b = [16.2, 17.3, 16.6]; % 成本系数b c = [1000, 970, 700]; % 成本系数c eg = [0.91, 0.95, 0.98]; % 碳排放系数 GT_P_max = [150, 80]; % 燃气轮机最大出力 E_bata = 0.81; % 碳捕集效率 P_yita = [1.05, 1.05, 1.05]; % 再生塔和压缩机系数 P_lamda = [0.269, 0.269, 0.269]; % 单位捕碳量能耗 P_G_max = [400, 455, 200]; % 火电机组最大出力 P_G_min = [200, 120, 100]; % 火电机组最小出力 R_u = 4 * [50, 50, 25]; % 机组爬坡率 % 碳捕集系统参数 P_D1 = 0.5 * [20, 20, 15]; % 碳捕集固定能耗 P_D = repmat(P_D1', 1, T); % 扩展至24小时 M_MEA = 61.08; % M_MEA摩尔质量 M_co2 = 44; % CO2摩尔质量 sita = 0.4; % 再生塔解析量 CR = 30; % 溶液浓度(%) rou_R = 1.01; % 溶液密度 V_CR = [2000, 2000, 2000]; % 溶液储液装置容量 K_R = 1.17; % 溶剂成本系数 fai = 1.5; % 溶剂损耗系数 % 分段线性化参数 h_gn = 5; % 上分段数 k_gn = 5; % 下分段数 gn = 5; % 火电机组分段数 % 计算分段边界 hl1 = (P_G_max - P_G_min) / h_gn; % 上分段步长 kl1 = (P_G_min - P_G_max) / k_gn; % 下分段步长 (注意负号) gl1 = (P_G_max - P_G_min) / gn; % 机组分段步长 %% 决策变量定义 % 电力相关 GT_P = sdpvar(2, T, 'full'); % 燃气轮机电出力 P_w = sdpvar(1, T, 'full'); % 风电出力 P_G = sdpvar(3, T, 'full'); % 火电机组总出力 P_J = sdpvar(3, T, 'full'); % 火电机组净出力 EB = sdpvar(2, T, 'full'); % 电锅炉出力 % 热力相关 GT_H = sdpvar(2, T, 'full'); % 燃气轮机热出力 EB_H = sdpvar(2, T, 'full'); % 电锅炉热出力 % 碳捕集系统 P_gas = sdpvar(2, T, 'full'); % 天然气需求 E_G = sdpvar(3, T, 'full'); % 总碳排放 E_total_co2 = sdpvar(3, T, 'full'); % 捕获的总碳排放 E_CG = sdpvar(3, T, 'full'); % 待捕集CO2量 P_B = sdpvar(3, T, 'full'); % 机组运行能耗 V_CA = sdpvar(3, T, 'full'); % 溶液体积 V_FY = sdpvar(3, T, 'full'); % 富液体积 V_PY = sdpvar(3, T, 'full'); % 贫液体积 % 需求响应 P_tran = sdpvar(1, T, 'full'); % 可转移电负荷 P_cut = sdpvar(1, T, 'full'); % 可削减电负荷 P_DE = sdpvar(1, T, 'full'); % 响应后电负荷 H_tran = sdpvar(1, T, 'full'); % 可转移热负荷 H_cut = sdpvar(1, T, 'full'); % 可削减热负荷 H_DE = sdpvar(1, T, 'full'); % 响应后热负荷 % 绝对值辅助变量 (用于线性化目标函数中的绝对值) P_tran_pos = sdpvar(1, T, 'full'); % P_tran正部 P_tran_neg = sdpvar(1, T, 'full'); % P_tran负部 H_tran_pos = sdpvar(1, T, 'full'); % H_tran正部 H_tran_neg = sdpvar(1, T, 'full'); % H_tran负部 % 分段线性化变量 w = binvar(gn, T, 3); % 分段激活标志 (机组×时间×分段) lambda = sdpvar(gn+1, T, 3);% 分段权重 (分段点×时间×机组) %% 约束构建 Cons = []; % ===== 基础约束 ===== % 风电出力约束 Cons = [Cons, 0 <= P_w <= P_wt_pre]; % 电锅炉约束 Cons = [Cons, sum(EB,1) + P_w <= P_wt_pre]; % 火电机组出力约束 for i = 1:3 Cons = [Cons, P_G_min(i) <= P_G(i,:) <= P_G_max(i)]; end % 绝对值分解约束 Cons = [Cons, P_tran == P_tran_pos - P_tran_neg]; Cons = [Cons, P_tran_pos >= 0, P_tran_neg >= 0]; Cons = [Cons, H_tran == H_tran_pos - H_tran_neg]; Cons = [Cons, H_tran_pos >= 0, H_tran_neg >= 0]; % 碳捕集系统静态约束 for i = 1:3 for t = 1:T % 总碳排放 Cons = [Cons, E_G(i,t) == eg(i) * P_G(i,t)]; % 捕集量约束 Cons = [Cons, 0 <= E_total_co2(i,t) <= P_yita(i) * E_bata * eg(i) * P_G_max(i)]; Cons = [Cons, E_total_co2(i,t) == E_CG(i,t) + 0.25 * E_bata * eg(i) * (P_G(i,t) - P_J(i,t))]; % 能耗平衡 Cons = [Cons, P_B(i,t) == P_lamda(i) * E_total_co2(i,t)]; Cons = [Cons, P_G(i,t) == P_J(i,t) + P_D(i,t) + P_B(i,t)]; % 净出力范围 P_J_min = P_G_min(i) - P_lamda(i) * P_yita(i) * E_bata * eg(i) * P_G_max(i) - P_D(i,t); P_J_max = P_G_max(i) - P_D(i,t); Cons = [Cons, P_J_min <= P_J(i,t) <= P_J_max]; end end % 旋转备用约束 (每小时) for t = 1:T Cons = [Cons, sum(R_u) >= 0.01 * P_DE(t)]; end % ===== 分段线性化约束 ===== for i = 1:3 % 对每台机组 seg_points = linspace(P_G_min(i), P_G_max(i), gn+1); % 分段点 for t = 1:T % 每个时间段 % 凸组合约束 Cons = [Cons, sum(lambda(:,t,i)) == 1, lambda(:,t,i) >= 0]; % 功率表达式 Cons = [Cons, P_G(i,t) == seg_points * lambda(:,t,i)]; % 分段激活约束 for seg = 1:gn Cons = [Cons, lambda(seg,t,i) <= w(seg,t,i)]; if seg < gn Cons = [Cons, lambda(seg+1,t,i) <= w(seg,t,i) + w(seg+1,t,i)]; end end Cons = [Cons, sum(w(:,t,i)) == 1]; % 仅激活一个分段 end end % ===== 机组运行约束 ===== % 爬坡约束 for i = 1:3 for t = 1:T-1 ramp_down = -R_u(i); ramp_up = R_u(i); Cons = [Cons, ramp_down <= P_G(i,t+1) - P_G(i,t) <= ramp_up]; end end % 燃气轮机约束 for j = 1:2 % 出力范围 Cons = [Cons, 0 <= GT_P(j,:) <= GT_P_max(j)]; % 天然气转换 Cons = [Cons, P_gas(j,:) == GT_P(j,:) / 0.98]; Cons = [Cons, 0 <= P_gas(j,:) <= 130]; % 热电关系 Cons = [Cons, GT_H(j,:) == 0.6 * (GT_P(j,:) / 0.4)]; Cons = [Cons, EB_H(j,:) == 0.98 * EB(j,:)]; end % ===== 碳捕集系统动态约束 ===== V_PY0 = [1000, 1000, 1000]; % 贫液初始值 V_FY0 = [1000, 1000, 1000]; % 富液初始值 for i = 1:3 % 初始条件 Cons = [Cons, V_FY(i,1) == V_FY0(i) - V_CA(i,1)]; Cons = [Cons, V_PY(i,1) == V_PY0(i) + V_CA(i,1)]; % 动态方程 for t = 2:T Cons = [Cons, V_FY(i,t) == V_FY(i,t-1) - V_CA(i,t)]; Cons = [Cons, V_PY(i,t) == V_PY(i,t-1) + V_CA(i,t)]; end % 溶液体积计算 for t = 1:T Cons = [Cons, V_CA(i,t) == (M_MEA * E_CG(i,t)) / (M_co2 * sita * CR * rou_R)]; Cons = [Cons, 0 <= V_FY(i,t) <= V_CR(i)]; Cons = [Cons, 0 <= V_PY(i,t) <= V_CR(i)]; end end % ===== 需求响应约束 ===== % 电负荷响应 for t = 1:T if ismember(t, [11,12,18,19,20]) % 高峰时段 Cons = [Cons, -0.02 * P_Load_e(t) <= P_cut(t) <= 0]; else Cons = [Cons, P_cut(t) == 0]; end end % 热负荷响应 for t = 1:T if ismember(t, [1,2,3,11,12]) % 高峰时段 Cons = [Cons, -0.02 * P_Load_h(t) <= H_cut(t) <= 0]; else Cons = [Cons, H_cut(t) == 0]; end end % 响应平衡约束 Cons = [Cons, -0.08 * P_Load_e <= P_tran <= 0.08 * P_Load_e]; Cons = [Cons, sum(P_tran) == 0]; Cons = [Cons, P_DE == P_Load_e + P_tran + P_cut]; Cons = [Cons, -0.08 * P_Load_h <= H_tran <= 0.08 * P_Load_h]; Cons = [Cons, sum(H_tran) == 0]; Cons = [Cons, H_DE == P_Load_h + H_tran + H_cut]; % ===== 系统平衡约束 ===== Cons = [Cons, P_w + sum(GT_P,1) + sum(P_J,1) == P_DE]; % 电力平衡 Cons = [Cons, sum(EB_H,1) + sum(GT_H,1) == H_DE]; % 热力平衡 %% 目标函数构建 % 需求响应成本 cost_DR_elec = 0.25 * gama_elec .* (P_tran_pos + P_tran_neg) + ... 0.65 * gama_elec .* (-P_cut); % P_cut<=0 cost_DR_heat = 0.25 * (H_tran_pos + H_tran_neg) + ... 0.65 * (-H_cut); % H_cut<=0 % 发电成本 cost_gen = zeros(1,T); for i = 1:3 cost_gen = cost_gen + a(i)*P_G(i,:).^2 + b(i)*P_G(i,:) + c(i); end % 碳捕集成本 cost_carbon = zeros(1,T); for i = 1:3 cost_carbon = cost_carbon + 1.5 * K_R * E_total_co2(i,:); end % 排放成本 cost_emission = zeros(1,T); for i = 1:3 net_emission = E_G(i,:) - E_total_co2(i,:); cost_emission = cost_emission + 12 * net_emission; end % 其他成本 cost_gas = 50 * sum(P_gas,1); cost_EB = gama_elec .* sum(EB,1); wind_curtail = max(0, P_wt_pre - P_w - sum(EB,1)); % 弃风量 cost_wind_curtail = 0.3 * gama_elec .* wind_curtail; % 总目标函数 F = sum(cost_DR_elec) + sum(cost_DR_heat) + sum(cost_gen) + ... sum(cost_carbon) + sum(cost_emission) + sum(cost_gas) + ... sum(cost_EB) + sum(cost_wind_curtail); %% 模型求解 ops = sdpsettings('solver', 'cplex', 'verbose', 2, 'usex0', 0); ops.cplex.mip.tolerances.mipgap = 1e-6; result = optimize(Cons, F, ops); %% 结果分析与可视化 if result.problem == 0 disp(['最优值: ', num2str(value(F))]); % 结果提取 val = @(x) value(x); % 简写取值函数 [GT_P_val, P_w_val, P_G_val] = deal(val(GT_P), val(P_w), val(P_G)); [EB_val, GT_H_val, EB_H_val] = deal(val(EB), val(GT_H), val(EB_H)); [P_J_val, E_G_val, E_total_co2_val] = deal(val(P_J), val(E_G), val(E_total_co2)); [P_DE_val, H_DE_val] = deal(val(P_DE), val(H_DE)); % 可视化设置 t = 1:T; colors = lines(7); % 获取7种区分度高的颜色 % 1. 系统功率平衡 figure('Name', '系统功率平衡', 'Position', [100, 100, 1200, 800]); % 电功率平衡 subplot(2,1,1); hold on; area(t, [val(P_J(1,:)); val(P_J(2,:)); val(P_J(3,:)); val(P_w); val(GT_P(1,:)); val(GT_P(2,:))]', ... 'LineWidth', 1.5); plot(t, P_Load_e, 'k--o', 'LineWidth', 2, 'MarkerSize', 8); plot(t, P_DE_val, 'm-*', 'LineWidth', 2, 'MarkerSize', 8); legend('火电1', '火电2', '火电3', '风电', '燃机1', '燃机2', ... '原始电负荷', '响应后电负荷', 'Location', 'best'); title('电功率平衡'); xlabel('时间 (h)'); ylabel('功率 (kW)'); grid on; box on; % 热功率平衡 subplot(2,1,2); hold on; area(t, [val(GT_H(1,:)); val(GT_H(2,:)); val(EB_H(1,:)); val(EB_H(2,:))]', ... 'LineWidth', 1.5); plot(t, P_Load_h, 'k--o', 'LineWidth', 2, 'MarkerSize', 8); plot(t, H_DE_val, 'm-*', 'LineWidth', 2, 'MarkerSize', 8); legend('燃机1热', '燃机2热', '电锅炉1', '电锅炉2', ... '原始热负荷', '响应后热负荷', 'Location', 'best'); title('热功率平衡'); xlabel('时间 (h)'); ylabel('功率 (kW)'); grid on; box on; % 2. 成本分解 figure('Name', '成本分解', 'Position', [100, 100, 1000, 600]); costs = [sum(val(cost_DR_elec)), sum(val(cost_DR_heat)), sum(val(cost_gen)), ... sum(val(cost_carbon)), sum(val(cost_emission)), ... sum(val(cost_gas)), sum(val(cost_EB)), sum(val(cost_wind_curtail))]; labels = {'需求响应(电)', '需求响应(热)', '发电成本', '碳捕集', ... '碳排放', '天然气', '电锅炉', '弃风惩罚'}; pie(costs, labels); title('总成本构成'); % 3. 碳捕集系统分析 figure('Name', '碳捕集分析', 'Position', [100, 100, 1200, 800]); % CO2捕集量 subplot(2,1,1); plot(t, val(E_total_co2(1,:)), 'r-o', 'LineWidth', 1.5); hold on; plot(t, val(E_total_co2(2,:)), 'b-s', 'LineWidth', 1.5); plot(t, val(E_total_co2(3,:)), 'g-^', 'LineWidth', 1.5); legend('机组1', '机组2', '机组3', 'Location', 'best'); title('CO2捕集量'); xlabel('时间 (h)'); ylabel('捕集量 (kg)'); grid on; box on; % 溶液体积变化 subplot(2,1,2); plot(t, val(V_FY(1,:)), 'r--o', 'LineWidth', 1.5); hold on; plot(t, val(V_PY(1,:)), 'r-s', 'LineWidth', 1.5); plot(t, val(V_FY(2,:)), 'b--^', 'LineWidth', 1.5); plot(t, val(V_PY(2,:)), 'b-d', 'LineWidth', 1.5); legend('机组1富液', '机组1贫液', '机组2富液', '机组2贫液', 'Location', 'best'); title('溶液体积动态变化'); xlabel('时间 (h)'); ylabel('体积 (m³)'); grid on; box on; else error('优化失败: %s', result.info); end
最新发布
11-18
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值