采用二次规划求解电力系统经济调度问题
电力系统经济调度是电力系统运行中的核心问题，其目标是在满足各种约束条件下，最小化发电总成本。
1. 描述
电力系统经济调度问题的数学模型如下：
目标函数：
\[\min \sum_{i=1}^{N} (a_i P_i^2 + b_i P_i + c_i)
\]
约束条件：
- 功率平衡约束：\(\sum_{i=1}^{N} P_i = P_D + P_L\)
- 发电机出力上下限：\(P_i^{\min} \leq P_i \leq P_i^{\max}, \quad i=1,2,\ldots,N\)
- 发电机爬坡率约束（可选）：\(|P_i(t) - P_i(t-1)| \leq R_i\)
其中：
- \(P_i\) 是第i台发电机的有功出力
- \(a_i, b_i, c_i\) 是第i台发电机的成本系数
- \(P_D\) 是系统总负荷
- \(P_L\) 是网络损耗（通常可以近似表示）
2. MATLAB实现
2.1 基本经济调度（不考虑网损和爬坡率）
function economic_dispatch_simple()
% 简单经济调度问题（不考虑网损和爬坡率约束）
% 发电机参数 [a, b, c, Pmin, Pmax]
generators = [
0.00375, 2.00, 0, 50, 200; % 发电机1
0.01750, 1.75, 0, 20, 100; % 发电机2
0.06250, 1.00, 0, 15, 50; % 发电机3
];
% 总负荷需求
PD = 180; % MW
% 提取参数
a = generators(:, 1);
b = generators(:, 2);
c = generators(:, 3);
Pmin = generators(:, 4);
Pmax = generators(:, 5);
n = size(generators, 1); % 发电机数量
% 构建二次规划问题
% 目标函数: min 0.5*x'*H*x + f'*x
H = 2 * diag(a); % 二次项系数
f = b; % 一次项系数
% 约束条件: Aeq*x = beq, lb <= x <= ub
Aeq = ones(1, n); % 功率平衡约束
beq = PD;
lb = Pmin; % 下限
ub = Pmax; % 上限
% 求解二次规划问题
options = optimoptions('quadprog', 'Display', 'iter', 'Algorithm', 'interior-point-convex');
[P, fval, exitflag] = quadprog(H, f, [], [], Aeq, beq, lb, ub, [], options);
% 显示结果
if exitflag > 0
fprintf('经济调度求解成功！\n');
fprintf('总成本: %.2f $\n', fval + sum(c));
fprintf('总发电量: %.2f MW\n', sum(P));
fprintf('总负荷: %.2f MW\n', PD);
fprintf('\n各发电机出力:\n');
for i = 1:n
fprintf('发电机%d: %.2f MW (成本: %.2f $/h)\n', i, P(i), a(i)*P(i)^2 + b(i)*P(i) + c(i));
end
% 绘制成本函数曲线和最优解
plot_cost_functions(generators, P);
else
fprintf('经济调度求解失败！\n');
end
end
function plot_cost_functions(generators, P_opt)
% 绘制发电机成本函数曲线和最优解
figure;
hold on;
colors = ['r', 'g', 'b', 'm', 'c'];
for i = 1:size(generators, 1)
a = generators(i, 1);
b = generators(i, 2);
c = generators(i, 3);
Pmin = generators(i, 4);
Pmax = generators(i, 5);
% 生成功率范围
P_range = linspace(Pmin, Pmax, 100);
% 计算成本
cost = a * P_range.^2 + b * P_range + c;
% 绘制成本曲线
plot(P_range, cost, [colors(i) '-'], 'LineWidth', 1.5, 'DisplayName', sprintf('发电机%d', i));
% 标记最优解
cost_opt = a * P_opt(i)^2 + b * P_opt(i) + c;
plot(P_opt(i), cost_opt, [colors(i) 'o'], 'MarkerSize', 8, 'MarkerFaceColor', colors(i));
end
xlabel('出力 (MW)');
ylabel('成本 ($/h)');
title('发电机成本函数');
legend('show');
grid on;
hold off;
end
2.2 考虑网损的经济调度
function economic_dispatch_with_losses()
% 考虑网损的经济调度问题
% 发电机参数 [a, b, c, Pmin, Pmax]
generators = [
0.00375, 2.00, 0, 50, 200; % 发电机1
0.01750, 1.75, 0, 20, 100; % 发电机2
0.06250, 1.00, 0, 15, 50; % 发电机3
];
% 总负荷需求
PD = 180; % MW
% B系数矩阵（用于计算网损）
B = [
0.0001, 0.0000, 0.0000;
0.0000, 0.0001, 0.0000;
0.0000, 0.0000, 0.0001;
];
B0 = [0.0000, 0.0000, 0.0000]';
B00 = 0;
% 提取参数
a = generators(:, 1);
b = generators(:, 2);
c = generators(:, 3);
Pmin = generators(:, 4);
Pmax = generators(:, 5);
n = size(generators, 1); % 发电机数量
% 构建二次规划问题
H = 2 * diag(a);
f = b;
lb = Pmin;
ub = Pmax;
% 使用迭代法求解考虑网损的经济调度
max_iter = 20;
tolerance = 1e-6;
lambda = 10; % 初始拉格朗日乘子
P = zeros(n, 1); % 初始发电计划
for iter = 1:max_iter
% 计算惩罚因子（考虑网损）
PL = P' * B * P + B0' * P + B00; % 网损
penalty = 2 * B * P + B0; % 网损对功率的导数
% 构建优化问题
% 目标函数: min 0.5*x'*H*x + f'*x
% 约束条件: sum(P_i) = PD + PL
% 使用近似方法处理网损约束
Aeq = ones(1, n);
beq = PD + PL;
% 求解二次规划
options = optimoptions('quadprog', 'Display', 'off', 'Algorithm', 'interior-point-convex');
[P_new, ~, exitflag] = quadprog(H, f, [], [], Aeq, beq, lb, ub, [], options);
if exitflag <= 0
error('二次规划求解失败！');
end
% 检查收敛
if max(abs(P_new - P)) < tolerance
P = P_new;
break;
end
P = P_new;
% 更新拉格朗日乘子（可选）
% 这里可以使用更复杂的方法更新lambda
fprintf('迭代 %d: 最大变化量 = %.6f\n', iter, max(abs(P_new - P)));
end
% 计算最终网损和总成本
PL = P' * B * P + B0' * P + B00;
total_cost = sum(a .* P.^2 + b .* P + c);
% 显示结果
fprintf('经济调度求解成功！\n');
fprintf('迭代次数: %d\n', iter);
fprintf('总成本: %.2f $\n', total_cost);
fprintf('总发电量: %.2f MW\n', sum(P));
fprintf('总负荷: %.2f MW\n', PD);
fprintf('网损: %.2f MW\n', PL);
fprintf('\n各发电机出力:\n');
for i = 1:n
fprintf('发电机%d: %.2f MW (成本: %.2f $/h)\n', i, P(i), a(i)*P(i)^2 + b(i)*P(i) + c(i));
end
% 绘制成本函数曲线和最优解
plot_cost_functions(generators, P);
end
2.3 考虑爬坡率约束的经济调度
function economic_dispatch_with_ramping()
% 考虑爬坡率约束的经济调度问题（多时段）
% 发电机参数 [a, b, c, Pmin, Pmax, RampUp, RampDown]
generators = [
0.00375, 2.00, 0, 50, 200, 40, 40; % 发电机1
0.01750, 1.75, 0, 20, 100, 30, 30; % 发电机2
0.06250, 1.00, 0, 15, 50, 20, 20; % 发电机3
];
% 负荷需求（多个时段）
PD = [180, 200, 170, 190, 210]; % MW
% 时段数量
T = length(PD);
% 提取参数
a = generators(:, 1);
b = generators(:, 2);
c = generators(:, 3);
Pmin = generators(:, 4);
Pmax = generators(:, 5);
RampUp = generators(:, 6);
RampDown = generators(:, 7);
n = size(generators, 1); % 发电机数量
% 构建多时段经济调度问题
% 决策变量: P(i,t) i=1..n, t=1..T
% 目标函数: min sum_{t=1}^T sum_{i=1}^n [a_i*P(i,t)^2 + b_i*P(i,t) + c_i]
H = 2 * kron(eye(T), diag(a)); % 二次项系数矩阵
f = repmat(b, T, 1); % 一次项系数向量
% 约束条件
% 1. 功率平衡约束: 每个时段 sum_i P(i,t) = PD(t)
Aeq = zeros(T, n*T);
for t = 1:T
Aeq(t, (t-1)*n+1:t*n) = ones(1, n);
end
beq = PD';
% 2. 发电机出力上下限
lb = repmat(Pmin, T, 1);
ub = repmat(Pmax, T, 1);
% 3. 爬坡率约束
A = [];
b_ramp = [];
for t = 2:T
for i = 1:n
% 上爬坡约束: P(i,t) - P(i,t-1) <= RampUp(i)
constraint = zeros(1, n*T);
constraint((t-1)*n+i) = 1;
constraint((t-2)*n+i) = -1;
A = [A; constraint];
b_ramp = [b_ramp; RampUp(i)];
% 下爬坡约束: P(i,t-1) - P(i,t) <= RampDown(i)
constraint = zeros(1, n*T);
constraint((t-1)*n+i) = -1;
constraint((t-2)*n+i) = 1;
A = [A; constraint];
b_ramp = [b_ramp; RampDown(i)];
end
end
% 合并所有不等式约束
A = [A; -eye(n*T)]; % 添加下界约束（转换为A*x <= b形式）
b = [b_ramp; -lb];
% 对于上界约束，我们需要单独处理
% 由于quadprog已经有ub参数，我们不需要在A中添加上界约束
% 求解二次规划问题
options = optimoptions('quadprog', 'Display', 'iter', 'Algorithm', 'interior-point-convex');
[P_all, fval, exitflag] = quadprog(H, f, A, b, Aeq, beq, [], ub, [], options);
% 重组结果
P = reshape(P_all, n, T);
% 显示结果
if exitflag > 0
fprintf('多时段经济调度求解成功！\n');
fprintf('总成本: %.2f $\n', fval + T*sum(c));
for t = 1:T
fprintf('\n时段 %d:\n', t);
fprintf('负荷需求: %.2f MW\n', PD(t));
fprintf('总发电量: %.2f MW\n', sum(P(:, t)));
for i = 1:n
fprintf('发电机%d: %.2f MW (成本: %.2f $/h)\n', i, P(i, t), ...
a(i)*P(i, t)^2 + b(i)*P(i, t) + c(i));
end
end
% 绘制多时段调度结果
plot_multiperiod_results(PD, P, generators);
else
fprintf('经济调度求解失败！\n');
end
end
function plot_multiperiod_results(PD, P, generators)
% 绘制多时段经济调度结果
[n, T] = size(P);
t = 1:T;
figure;
% 绘制各发电机出力和总负荷
subplot(2, 1, 1);
hold on;
% 计算累积出力
cumulative = zeros(T, n);
for i = 1:n
if i == 1
cumulative(:, i) = P(i, :)';
else
cumulative(:, i) = cumulative(:, i-1) + P(i, :)';
end
end
% 绘制堆叠区域图
colors = ['r', 'g', 'b', 'm', 'c'];
for i = n:-1:1
if i == 1
area(t, cumulative(:, i), 'FaceColor', colors(i), 'DisplayName', sprintf('发电机%d', i));
else
area(t, cumulative(:, i), 'FaceColor', colors(i), 'DisplayName', sprintf('发电机%d', i));
end
end
% 绘制总负荷
plot(t, PD, 'k--', 'LineWidth', 2, 'DisplayName', '负荷需求');
xlabel('时段');
ylabel('出力 (MW)');
title('多时段经济调度结果');
legend('show');
grid on;
hold off;
% 绘制各发电机成本
subplot(2, 1, 2);
hold on;
a = generators(:, 1);
b = generators(:, 2);
c = generators(:, 3);
cost = zeros(n, T);
for i = 1:n
cost(i, :) = a(i) * P(i, :).^2 + b(i) * P(i, :) + c(i);
plot(t, cost(i, :), [colors(i) '-o'], 'LineWidth', 1.5, 'DisplayName', sprintf('发电机%d', i));
end
xlabel('时段');
ylabel('成本 ($/h)');
title('各发电机成本变化');
legend('show');
grid on;
hold off;
end
2.4 综合经济调度（考虑阀点效应）
function economic_dispatch_valve_point()
% 考虑阀点效应的经济调度问题
% 发电机参数 [a, b, c, Pmin, Pmax, e, f] (e,f是阀点效应系数)
generators = [
0.00375, 2.00, 0, 50, 200, 0.063, 0.5; % 发电机1
0.01750, 1.75, 0, 20, 100, 0.098, 0.6; % 发电机2
0.06250, 1.00, 0, 15, 50, 0.120, 0.7; % 发电机3
];
% 总负荷需求
PD = 180; % MW
% 提取参数
a = generators(:, 1);
b = generators(:, 2);
c = generators(:, 3);
Pmin = generators(:, 4);
Pmax = generators(:, 5);
e = generators(:, 6);
f = generators(:, 7);
n = size(generators, 1); % 发电机数量
% 阀点效应成本函数: F_i(P_i) = a_i*P_i^2 + b_i*P_i + c_i + |e_i*sin(f_i*(P_i_min - P_i))|
% 这是一个非凸函数，不能直接用二次规划求解
% 我们可以使用分段线性近似或启发式算法
% 这里使用分段线性近似方法
% 为每台发电机创建分段点
num_segments = 20; % 每台发电机的分段数
P_segments = cell(n, 1);
cost_segments = cell(n, 1);
for i = 1:n
P_range = linspace(Pmin(i), Pmax(i), num_segments);
cost_valve = a(i)*P_range.^2 + b(i)*P_range + c(i) + ...
e(i)*abs(sin(f(i)*(Pmin(i) - P_range)));
P_segments{i} = P_range;
cost_segments{i} = cost_valve;
end
% 使用混合整数规划或线性规划求解
% 这里简化处理，使用fmincon求解
% 初始猜测（均匀分配）
P0 = ones(n, 1) * PD / n;
% 定义目标函数
objective = @(P) sum(a .* P.^2 + b .* P + c + e .* abs(sin(f .* (Pmin - P))));
% 约束条件
Aeq = ones(1, n);
beq = PD;
lb = Pmin;
ub = Pmax;
% 求解
options = optimoptions('fmincon', 'Display', 'iter', 'Algorithm', 'interior-point');
[P, fval, exitflag] = fmincon(objective, P0, [], [], Aeq, beq, lb, ub, [], options);
% 显示结果
if exitflag > 0
fprintf('考虑阀点效应的经济调度求解成功！\n');
fprintf('总成本: %.2f $\n', fval);
fprintf('总发电量: %.2f MW\n', sum(P));
fprintf('总负荷: %.2f MW\n', PD);
fprintf('\n各发电机出力:\n');
for i = 1:n
cost_i = a(i)*P(i)^2 + b(i)*P(i) + c(i) + e(i)*abs(sin(f(i)*(Pmin(i) - P(i))));
fprintf('发电机%d: %.2f MW (成本: %.2f $/h)\n', i, P(i), cost_i);
end
% 绘制成本函数曲线和最优解
plot_valve_point_functions(generators, P);
else
fprintf('经济调度求解失败！\n');
end
end
function plot_valve_point_functions(generators, P_opt)
% 绘制考虑阀点效应的发电机成本函数曲线和最优解
figure;
hold on;
colors = ['r', 'g', 'b', 'm', 'c'];
for i = 1:size(generators, 1)
a = generators(i, 1);
b = generators(i, 2);
c = generators(i, 3);
Pmin = generators(i, 4);
Pmax = generators(i, 5);
e = generators(i, 6);
f = generators(i, 7);
% 生成功率范围
P_range = linspace(Pmin, Pmax, 1000);
% 计算成本（考虑阀点效应）
cost = a * P_range.^2 + b * P_range + c + e * abs(sin(f * (Pmin - P_range)));
% 绘制成本曲线
plot(P_range, cost, [colors(i) '-'], 'LineWidth', 1.5, 'DisplayName', sprintf('发电机%d', i));
% 标记最优解
cost_opt = a * P_opt(i)^2 + b * P_opt(i) + c + e * abs(sin(f * (Pmin - P_opt(i))));
plot(P_opt(i), cost_opt, [colors(i) 'o'], 'MarkerSize', 8, 'MarkerFaceColor', colors(i));
end
xlabel('出力 (MW)');
ylabel('成本 ($/h)');
title('考虑阀点效应的发电机成本函数');
legend('show');
grid on;
hold off;
end
3. 主函数和示例运行
function main_economic_dispatch()
% 主函数：运行各种经济调度示例
fprintf('=== 电力系统经济调度问题求解 ===\n\n');
% 简单经济调度
fprintf('1. 简单经济调度（不考虑网损和爬坡率）\n');
economic_dispatch_simple();
fprintf('\n\n');
% 考虑网损的经济调度
fprintf('2. 考虑网损的经济调度\n');
economic_dispatch_with_losses();
fprintf('\n\n');
% 考虑爬坡率约束的经济调度
fprintf('3. 考虑爬坡率约束的多时段经济调度\n');
economic_dispatch_with_ramping();
fprintf('\n\n');
% 考虑阀点效应的经济调度
fprintf('4. 考虑阀点效应的经济调度\n');
economic_dispatch_valve_point();
fprintf('\n=== 所有示例运行完成 ===\n');
end
4. 如何运行
要运行这些经济调度示例，只需在MATLAB命令窗口中调用主函数：
main_economic_dispatch();
或者直接运行特定的经济调度函数：
% 运行简单经济调度
economic_dispatch_simple();
% 运行考虑网损的经济调度
economic_dispatch_with_losses();
% 运行考虑爬坡率约束的经济调度
economic_dispatch_with_ramping();
% 运行考虑阀点效应的经济调度
economic_dispatch_valve_point();
推荐代码 采用二次规划求解电力系统经济调度问题 www.youwenfan.com/contentcng/50899.html
5. 结果分析
运行上述代码后，你将得到：
- 简单经济调度：各发电机的最优出力分配和总成本，以及成本函数曲线图。
- 考虑网损的经济调度：考虑网络损耗后的最优发电计划，显示网损值和迭代过程。
- 考虑爬坡率约束的多时段经济调度：多个时段的发电计划，显示各时段各发电机的出力和成本变化。
- 考虑阀点效应的经济调度：考虑阀点效应后的最优发电计划，显示非凸成本函数曲线。
6. 扩展功能
你可以进一步扩展这个实现，例如：
- 添加更多约束：如输电线路容量约束、备用容量约束等。
- 考虑可再生能源：如风能、太阳能的随机性和不确定性。
- 动态经济调度：考虑机组启停成本和最小启停时间约束。
- 市场环境下的经济调度：考虑电力市场竞价机制。
- 多目标优化：同时考虑经济性和环保性（减排）。
这个实现提供了一个完整的电力系统经济调度求解框架，你可以根据需要修改发电机参数、负荷数据和约束条件。
浙公网安备 33010602011771号