P1734_最大约数和

题目描述

选取和不超过S的若干个不同的正整数,使得所有数的约数(不含它本身)之和最大。

输入输出格式
输入格式:
输入一个正整数S。

输出格式:
输出最大的约数之和。

输入输出样例
输入样例#1:11
输出样例#1:9

说明

  • 样例说明 取数字4和6,可以得到最大值(1+2)+(1+2+3)=9。
  • 数据规模S<=1000

    和不超过某个数,又要让约数和最大,这很容易让人想起背包问题。想到这里问题就解决了。
    把每个数的大小看作weight数组,每个数本身的约数和看作value数组。这就是个01背包问题。

#include <cstdio>
#include <iostream>
#include <algorithm>

using namespace std;

int s, v[1001], dp[1001];

int num(
检查优化重构以下代码:%% 日前调度优化模型 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('体积 ()'); grid on; box on; else error('优化失败: %s', result.info); end
最新发布
11-18
再次优化以下代码,重构后输出:%% 日前调度优化模型 clc; clear; %% 参初始化 T = 24; % 调度周期 % 预测据 P_prew = 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]; EleLoad = 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]; HeatLoad = 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; % 火电机组分段 % 初始化分段边界 [h_min, h_max] = deal([200, 120, 100], [401, 456, 201]); [k_min, k_max] = deal([-400, -455, -200], [-199, -119, -99]); % 计算分段边界 hl1 = (h_max - h_min) / h_gn; kl1 = (k_max - k_min) / k_gn; gl1 = (P_G_max - P_G_min) / gn; hl2 = arrayfun(@(min, max, step) min:step:max, h_min, h_max, hl1, 'UniformOutput', false); kl2 = arrayfun(@(min, max, step) min:step:max, k_min, k_max, kl1, 'UniformOutput', false); gl2 = arrayfun(@(min, max, step) min:step:max, P_G_min, P_G_max, gl1, 'UniformOutput', false); %% 决策变量定义 % 电力相关 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'); % 响应后热负荷 % 分段线性化变量 % 火电机组分段 w = sdpvar(gn+1, T, 3); % 三维组: [分段点, 时间, 机组] z = binvar(gn, T, 3); % 三维二进制变量 P_G_line = sdpvar(3, T, 'full'); % 上/下分段变量 h_w = sdpvar(h_gn+1, T, 3); h_z = binvar(h_gn, T, 3); k_w = sdpvar(k_gn+1, T, 3); k_z = binvar(k_gn, T, 3); y1 = sdpvar(3, T, 'full'); y2 = sdpvar(3, T, 'full'); %% 约束构建 C = []; % ===== 基础约束 ===== for t = 1:T % 风电出力约束 C = [C, 0 <= P_w(t) <= P_prew(t)]; % 电锅炉约束 C = [C, sum(EB(:, t)) + P_w(t) <= P_prew(t)]; for i = 1:3 % 火电机组出力约束 C = [C, P_G_min(i) <= P_G(i, t) <= P_G_max(i)]; % 碳捕集系统约束 C = [C, E_G(i, t) == eg(i) * P_G(i, t)]; % 总碳排放 C = [C, E_total_co2(i, t) == E_CG(i, t) + 0.25 * E_bata * eg(i) * (y1(i, t) - y2(i, t))]; % 捕获总量 C = [C, 0 <= E_total_co2(i, t) <= P_yita(i) * E_bata * eg(i) * P_G_max(i)]; C = [C, P_B(i, t) == P_lamda(i) * E_total_co2(i, t)]; % 运行能耗 C = [C, P_G(i, t) == P_J(i, t) + P_D(i, t) + P_B(i, t)]; % 总功率平衡 C = [C, P_G_min(i) - P_lamda(i) * P_yita(i) * E_bata * eg(i) * P_G_max(i) - P_D(i, t) <= P_J(i, t) <= ... P_G_max(i) - P_D(i, t)]; % 净出力范围 end end % 旋转备用约束 C = [C, min(sum(R_u), sum(P_G_max) - sum(P_G)) >= 0.01 * max(P_DE)]; % ===== 分段线性化约束 ===== % 通用分段约束函 function C = add_segment_constraints(C, w, z, values, var, T, n_segments) n_units = size(var, 1); for i = 1:n_units C = [C, var(i, :) == values{i} * squeeze(w(:, :, i))]; for t = 1:T % 凸组合约束 C = [C, sum(w(:, t, i)) == 1, w(:, t, i) >= 0]; % 相邻约束 C = [C, w(1, t, i) <= z(1, t, i)]; for seg = 2:n_segments C = [C, w(seg, t, i) <= z(seg-1, t, i) + z(seg, t, i)]; end C = [C, w(n_segments+1, t, i) <= z(n_segments, t, i)]; % 激活一个分段 C = [C, sum(z(:, t, i)) == 1]; end end end % 应用分段约束 C = add_segment_constraints(C, w, z, gl2, P_G_line, T, gn); C = add_segment_constraints(C, h_w, h_z, hl2, y1, T, h_gn); C = add_segment_constraints(C, k_w, k_z, kl2, y2, T, k_gn); % 连接变量 for i = 1:3 C = [C, P_G(i, :) == gl2{i} * squeeze(w(:, :, i))]; C = [C, P_G_line(i, :) == (gl2{i}.^2) * squeeze(w(:, :, i))]; end C = [C, y1 >= y2]; % 上下界关系约束 % ===== 机组运行约束 ===== % 爬坡约束 for t = 1:T-1 for i = 1:3 C = [C, -R_u(i) <= P_G(i, t+1) - P_G(i, t) <= R_u(i)]; end for j = 1:2 C = [C, -50 <= GT_P(j, t+1) - GT_P(j, t) <= 50]; C = [C, -30 <= EB(j, t+1) - EB(j, t) <= 30]; end end % 燃气轮机约束 for t = 1:T for j = 1:2 C = [C, 0 <= P_gas(j, t) <= 130]; C = [C, 0 <= GT_P(j, t) <= GT_P_max(j)]; C = [C, GT_P(j, t) == 0.98 * P_gas(j, t)]; C = [C, GT_H(j, t) == 0.6 * (GT_P(j, t) / 0.4)]; C = [C, 0 <= EB(j, t) <= P_prew(t) - P_w(t)]; C = [C, EB_H(j, t) == 0.98 * EB(j, t)]; end end % ===== 碳捕集系统动态约束 ===== V_PY0 = [1000, 1000, 1000]; % 贫液初始值 V_FY0 = [1000, 1000, 1000]; % 富液初始值 for i = 1:3 C = [C, V_PY0(i) == V_PY(i, T)]; % 终值平衡 C = [C, V_FY0(i) == V_FY(i, T)]; % 终值平衡 for t = 1:T C = [C, V_CA(i, t) == (M_MEA * E_CG(i, t)) / (M_co2 * sita * CR * rou_R)]; % 溶液体积 C = [C, 0 <= V_FY(i, t) <= V_CR(i)]; C = [C, 0 <= V_PY(i, t) <= V_CR(i)]; C = [C, 0 <= (0.25 * eg(i) * (y1(i, t) - y2(i, t))) <= P_G(i, t)]; end for t = 1:T-1 C = [C, V_FY(i, t+1) == V_FY(i, t) - V_CA(i, t+1)]; % 富液动态 C = [C, V_PY(i, t+1) == V_PY(i, t) + V_CA(i, t+1)]; % 贫液动态 end end % ===== 需求响应约束 ===== % 电负荷响应 for t = 1:T if (t >= 11 && t <= 12) || (t >= 18 && t <= 20) C = [C, -0.02 * EleLoad(t) <= P_cut(t) <= 0]; else C = [C, P_cut(t) == 0]; end end % 热负荷响应 for t = 1:T if (t >= 0 && t <= 3) || (t >= 11 && t <= 12) C = [C, -0.02 * HeatLoad(t) <= H_cut(t) <= 0]; else C = [C, H_cut(t) == 0]; end end % 响应平衡约束 C = [C, -0.08 * EleLoad <= P_tran <= 0.08 * EleLoad]; C = [C, sum(P_tran) == 0]; C = [C, P_DE == EleLoad + P_tran + P_cut]; C = [C, -0.08 * HeatLoad <= H_tran <= 0.08 * HeatLoad]; C = [C, sum(H_tran) == 0]; C = [C, H_DE == HeatLoad + H_tran + H_cut]; % ===== 系统平衡约束 ===== C = [C, P_w + sum(GT_P, 1) + sum(P_J, 1) == P_DE]; % 电力平衡 C = [C, sum(EB_H, 1) + sum(GT_H, 1) == H_DE]; % 热力平衡 %% 目标函构建 % 组件成本计算 cost_DR_elec = 0.25 * gama_elec .* abs(P_tran) + 0.65 * gama_elec .* abs(P_cut); cost_DR_heat = 0.25 * abs(H_tran) + 0.65 * abs(H_cut); cost_gen = zeros(1, T); cost_carbon = zeros(1, T); cost_emission = zeros(1, T); for i = 1:3 cost_gen = cost_gen + a(i) * P_G_line(i, :) + b(i) * P_G(i, :) + c(i); cost_carbon = cost_carbon + 1.5 * K_R * (0.25 * E_bata * eg(i) * (y1(i, :) - y2(i, :)) + E_CG(i, :)); cost_emission = cost_emission + 12 * (E_G(i, :) - 0.25 * E_bata * eg(i) * (y1(i, :) - y2(i, :)) - E_CG(i, :) - 0.7 * P_G(i, :)); end cost_gas = 50 * sum(P_gas, 1); cost_EB = gama_elec .* sum(EB, 1); cost_wind_curtail = 0.3 * gama_elec .* (P_prew - P_w - sum(EB, 1)); % 总目标函 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(C, F, ops); if result.problem == 0 disp(['最优值: ', num2str(value(F))]); % 结果提取 [GT_P_val, P_w_val, P_G_val, EB_val] = deal(value(GT_P), value(P_w), value(P_G), value(EB)); [GT_H_val, EB_H_val, P_gas_val] = deal(value(GT_H), value(EB_H), value(P_gas)); [E_G_val, E_total_co2_val, E_CG_val] = deal(value(E_G), value(E_total_co2), value(E_CG)); [P_B_val, P_J_val, V_CA_val] = deal(value(P_B), value(P_J), value(V_CA)); [y1_val, y2_val] = deal(value(y1), value(y2)); [P_tran_val, P_cut_val, P_DE_val] = deal(value(P_tran), value(P_cut), value(P_DE)); [H_tran_val, H_cut_val, H_DE_val] = deal(value(H_tran), value(H_cut), value(H_DE)); %% 结果可视化 % 原始据图 figure('Name', '原始据'); plot(1:T, P_prew, 'r-o', 'LineWidth', 1.5, 'DisplayName', '风电出力'); hold on; plot(1:T, EleLoad, 'g-*', 'LineWidth', 1.5, 'DisplayName', '电负荷'); plot(1:T, HeatLoad, 'b-^', 'LineWidth', 1.5, 'DisplayName', '热负荷'); legend show; xlabel('时间/h'); ylabel('功率/kW'); title('原始据'); grid on; % 电力平衡图 Plot_EleNet = [P_J_val(1, :); P_J_val(2, :); P_J_val(3, :); P_w_val + sum(EB_val, 1); GT_P_val(1, :); GT_P_val(2, :); sum(EB_val, 1)]'; figure('Name', '电功率平衡'); bar(1:T, Plot_EleNet, 'stacked'); hold on; plot(1:T, EleLoad, 'k-', 'LineWidth', 2, 'DisplayName', '原始电负荷'); plot(1:T, P_DE_val, 'r-o', 'LineWidth', 1.5, 'DisplayName', '响应后电负荷'); legend([{'火电1','火电2','火电3','风电','燃气轮机1','燃气轮机2','电锅炉'}, {'原始电负荷','响应后电负荷'}]); xlabel('时间/h'); ylabel('功率/kW'); title('电功率平衡'); grid on; % 热力平衡图 Plot_HeatNet = [sum(GT_H_val, 1); sum(EB_H_val, 1)]'; figure('Name', '热功率平衡'); bar(1:T, Plot_HeatNet, 'stacked'); hold on; plot(1:T, HeatLoad, 'k-', 'LineWidth', 2, 'DisplayName', '原始热负荷'); plot(1:T, H_DE_val, 'r-o', 'LineWidth', 1.5, 'DisplayName', '响应后热负荷'); legend('燃气轮机', '电锅炉', '原始热负荷', '响应后热负荷'); xlabel('时间/h'); ylabel('功率/kW'); title('热功率平衡'); grid on; % 风电出力图 figure('Name', '风电功率'); plot(1:T, P_w_val + sum(EB_val, 1), 'g-o', 'LineWidth', 1.5, 'DisplayName', '实际利用'); hold on; plot(1:T, P_prew, 'k-', 'LineWidth', 1.5, 'DisplayName', '预测出力'); legend show; xlabel('时间/h'); ylabel('功率/kW'); title('风电功率'); grid on; % 碳捕集能耗图 figure('Name', '碳捕集能耗'); plot(1:T, sum(P_B_val, 1), 'g-^', 'LineWidth', 1.5); xlabel('时间/h'); ylabel('功率/kW'); title('碳捕集能耗'); grid on; % 火电机组净出力 figure('Name', '火电机组净出力'); plot(1:T, sum(P_J_val, 1), 'r-o', 'LineWidth', 1.5); xlabel('时间/h'); ylabel('功率/kW'); title('火电机组净出力'); grid on; % 二氧化碳捕集量 figure('Name', '碳捕集量'); plot(1:T, E_total_co2_val(1, :), 'r-', 'LineWidth', 1.5, 'DisplayName', '机组1'); hold on; plot(1:T, E_total_co2_val(2, :), 'g-^', 'LineWidth', 1.5, 'DisplayName', '机组2'); plot(1:T, E_total_co2_val(3, :), 'b-*', 'LineWidth', 1.5, 'DisplayName', '机组3'); legend show; xlabel('时间/h'); ylabel('碳捕集量'); title('机组二氧化碳捕集量'); grid on; else error('优化失败: %s', result.info); end
11-18
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值