电力系统优化技术
Economic Dispatch(5.7 class)
问题
问题一
拉格朗日变换和线性规划的关系
这里讲到的经济调度（ED）主要使用到的方法为拉格朗日变换，无非就是给P加上各种限制（如$\sum_{i=1}^{NG}P_{ai}=P_{D}$、$P_{min,i}$ $<=$ $P_ai$ ）等，这其实跟运筹学里的线性规划很类似，但是我翻看之前线性规划的笔记却没有查看到拉格朗日变换
线性规划为什么不用拉格朗日？ 线性规划有不等式约束 $x \ge 0$，而基础拉格朗日法只处理等式约束。
而实际上，线性规划的对偶变量（影子价格）和经济调度的 $\lambda$，本质上是同一个东西
问题二
最开头的F（P）的模型怎么得到的
另外，我真的很好奇这个$F(P_{ai})$的模型是怎么建出来的，这个$P_{ai}$真的只是我们通常认为的成本吗？或者说这个公式的实际意义是什么？
$F(P) = \rho \cdot H(P)$ 的物理含义：
- $H(P)$ — 燃料消耗函数，描述的是：发电机输出 $P$ (MW) 电力时，每小时需要烧多少燃料。来自厂家提供的热耗率曲线，通过实际测试得到
- 乘以燃料价格 $\rho$， $\rho$ 是燃料单价，由市场决定
问题三
为什么能够使用拉格朗日变换来解决这类问题
经济调度问题是一个等式约束下的最优化问题，而拉格朗日乘数法正是处理这类问题的标准数学工具
这里强调下拉格朗日乘数法的核心思想：把约束条件”吞进”目标函数里，把有约束问题变成无约束问题
补充下：拉格朗日乘子 $\lambda$ 在这个问题里有真实的物理含义， 就是系统增量（边际）成本，即多发 1MW 负荷需要多花多少钱。
问题四
IC1、IC2到底是什么？等微增原则要怎么理解，尤其是“所有单元有一样的增量成本”里面的“单元”指的是什么？
实际上：IC即为单台电机的增量成本（或者说边际成本），通俗理解就是在当前出力水平上再多发 1MW，每小时要多花多少钱。
而$\lambda$是系统增量成本。理论上，最优解情况下全部的IC都要等于$\lambda$，这就得到了等微增原则。
首先，等微增原则中的“单元”就是单台发电机组，当全部IC 全等于 $\lambda，任何微调都无法再降低成本，这就是最优解
问题五
考虑系统损耗的情况：这里PL的系数B到底是啥？矩阵？
问题六
最后讲的近似计算要怎么计算
见例题一
受限优化和最优化条件——拉格朗日数乘法
Economic Dispatch(经济调度)
最常见的模型（最简单）：
衍生出来的概念
- 增量成本IC/边际成本MC
- 系统增量成本
- 等微增原则
- 沉默成本
标准解法
进阶一：考虑发电容量约束
问题模型如下：
得到的拉格朗日变换：
解题步骤
要求$\lambda$和$P_D$的关系曲线
进阶二：考虑系统损耗（难）
其中的$$P_{L}(P_{G})$$
$$P_{L}=B_{0}+B_{1}^{T}P_{G}+\frac{1}{2}P_{G}^{T}B_{2}P_{G}$$
Example
$$P_{G}=\begin{bmatrix}P_{G1}\\P_{G2}\\P_{G3}\end{bmatrix}$$
$$B_{1}=\begin{bmatrix}B_{1,1}\\B_{1,2}\\B_{1,3}\end{bmatrix}$$
$$B_{2}=\begin{bmatrix}B_{2,11}&B_{2,12}&B_{2,13}\\B_{2,21}&B_{2,22}&B_{2,23}\\B_{2,31}&B_{2,32}&B_{2,33}\end{bmatrix}$$
$$\text{B-coefficient method}$$
拉格朗日变换以及重要参数
理想解题步骤
实际方法
例题一
解答
首先先用二的表格，计算忽略损耗的PG
再将得到的PG分别填入表格五（系统损耗），不断迭代————填入的a、b、c、PD、c不变，不断将产生的PGi新填回到黄色的PGi中，直到差值很小时
Linear Programming（线性规划）5.12
从经济调度到线性规划
上面我们讨论的经济调度（ED），核心方法是拉格朗日乘子法——把等式约束”吞进”目标函数，转化为无约束优化。但我们也提到了一个问题：线性规划为什么不用拉格朗日？
答案其实很简单：拉格朗日乘子法基础形式只管等式约束，而线性规划的”变量 ≥ 0”是不等式约束。当然拉格朗日可以推广到 KKT 条件来处理不等式，但那已经是非线性规划的范畴了。
不过，两者在本质上是相通的——线性规划的对偶变量（也叫影子价格，shadow price）和经济调度中的 $\lambda$，是同一个东西：都表示”放宽一个单位的约束，目标函数能改善多少”。在电力系统中，这就是多发 1MW 需要多花多少钱。
那什么时候我们需要线性规划？很简单：当成本函数是线性的、或者我们可以用分段线性来近似时，经济调度本身就是线性规划问题。
这篇笔记的理论部分主要来自运筹学课程（梁军老师），但我们会在最后回到电力系统的应用中去。
什么是线性规划
线性规划（Linear Programming, LP）研究的是：在一组线性约束条件下，如何决策使得一个线性目标函数达到最优（最大或最小）。
三个基本要素：
- 决策变量：我们要决定的量，比如发电机的出力 $P_1, P_2, \dots$
- 目标函数：我们想最大化或最小化的东西，比如总成本 $c_1 P_1 + c_2 P_2$
- 约束条件：决策变量必须满足的限制，比如发电容量上下限、功率平衡
电力系统中的简单例子：两台发电机，成本系数分别为 $c_1=2$ 和 $c_2=5$（元/MWh），有功出力 $P_1, P_2$（MW），总负荷 100MW，容量约束分别为 $[0,80]$ 和 $[0,60]$。目标是总成本最低：
$$\min z = 2P_1 + 5P_2$$
$$s.t. \quad P_1 + P_2 = 100, \quad 0 \le P_1 \le 80, \quad 0 \le P_2 \le 60$$
这就构成了一个线性规划问题。
经典例子（来自 PPT）：某工厂计划生产甲、乙两种产品。已知生产单位产品所需的电力消耗、污染指数的限制以及利润如下表：
| 甲 | 乙 | 相关限制 | |
|---|---|---|---|
| 耗电（千瓦） | 1 | 1 | 6 |
| 污染指数 | -1 | 2 | 8 |
| 利润（万元） | 3 | 1 |
此外，乙产品的产量必须是甲产品的 2 倍以上。设 $x_1, x_2$ 为甲、乙产量，$z$ 为总利润：
$$\max z = 3x_1 + x_2$$
$$s.t. \quad
\begin{cases}
x_1 + x_2 \le 6 \cr
-x_1 + 2x_2 \le 8 \cr
2x_1 - x_2 \le 0 \cr
x_1, x_2 \ge 0
\end{cases}$$
图解法：两个变量的直观理解
当只有两个决策变量时，我们可以在平面上直接画出可行域和目标函数，直观理解线性规划的性质。
以上面的工厂问题为例：以 $x_1, x_2$ 为坐标轴，画出四条约束直线，它们围成的区域就是可行域（所有满足约束的解的集合）：
然后画目标函数的等值线 $3x_1 + x_2 = z$，沿 $z$ 增大的方向平移，当等值线离开可行域前碰到的最后一个点，就是最优解。
图解法揭示了四种可能的解：
唯一最优解
目标函数等值线与可行域恰好在一个顶点处”擦肩而过”。这是最常见的情况。比如上面的工厂问题，最优解在可行域的某个顶点上。无穷多最优解
当目标函数等值线的斜率与某条约束边界平行时，整条边上的所有点都是最优解。比如把目标函数改成 $\max z = x_1 + x_2$，等值线 $x_1 + x_2 = z$ 与约束 $x_1 + x_2 \le 6$ 的边界平行，那么这条边上所有点都是最优解。无界解
可行域在目标函数增大的方向上无限延伸（无界可行域），目标值可以无限变大。通常是遗漏了某个关键约束条件。例如如果去掉 $x_1 + x_2 \le 6$ 这条约束，可行域向右上方无限延伸，利润就可以无限增大了。无可行解
约束条件之间存在矛盾，可行域为空集。比如要求 $x_1 + x_2 \ge 6$ 且 $2x_1 + 2x_2 \le 8$，显然两个半平面没有交集。
图解法的启示：
- 可行域是一个凸集（凸多边形）
- 最优解若存在，一定能在可行域的顶点上找到
- 但对于三维及以上的问题，我们需要代数方法——这就是单纯形法
标准形式与转换规则
为了用代数方法统一求解，我们需要把所有线性规划问题转化为标准形。
标准形：
$$\max z = CX$$
$$s.t. \quad AX = b, \quad X \ge 0$$
三个关键词：目标函数是 max、约束是等式、变量都 ≥ 0。
从一般形式转化为标准形，需要掌握以下规则：
| 情况 | 转化方法 |
|---|---|
| $\min z$ | 令 $z’ = -z$，改为 $\max z’$ |
| $\sum a_{ij}x_j \le b_i$ | 加松弛变量 $x_{si} \ge 0$，变为 $\sum a_{ij}x_j + x_{si} = b_i$ |
| $\sum a_{ij}x_j \ge b_i$ | 减剩余变量 $x_{si} \ge 0$，变为 $\sum a_{ij}x_j - x_{si} = b_i$ |
| $x_j \le 0$ | 令 $x_j’ = -x_j \ge 0$ |
| $x_j$ 无约束（free） | 令 $x_j = x_j’ - x_j’’$，其中 $x_j’, x_j’’ \ge 0$ |
松弛变量和剩余变量在目标函数中的系数为 0（不影响原问题）。
转化例题：
将 $\min z = x_1 + 2x_2 - 3x_3$ 化为标准形：
$$s.t. \quad
\begin{cases}
x_1 + x_2 + x_3 \le 7 \cr
x_1 + x_2 + x_3 \ge 2 \cr
-3x_1 + x_2 + 2x_3 = 5 \cr
x_1, x_2 \ge 0, x_3 \text{ 无约束}
\end{cases}$$
解：
- $x_3$ 无约束：令 $x_3 = x_4 - x_5$, $x_4, x_5 \ge 0$
- $\le 7$：加松弛变量 $x_6 \ge 0$：$x_1 + x_2 + (x_4 - x_5) + x_6 = 7$
- $\ge 2$：减剩余变量 $x_7 \ge 0$：$x_1 + x_2 + (x_4 - x_5) - x_7 = 2$
- 目标函数：$\min z \rightarrow \max z’ = -z = x_1 - 2x_2 + 3x_4 - 3x_5 + 0x_6 + 0x_7$
注意：松弛/剩余变量在目标函数中系数为 0。
代数基础：基可行解与顶点
把线性规划的标准形写开：
$$\max z = CX, \quad s.t. \quad AX = b, \quad X \ge 0$$
其中 $A$ 是 $m \times n$ 矩阵，且 $m < n$（约束个数 < 变量个数，这才有多解的空间）。设 $A$ 的秩为 $m$。
基：从 $A$ 的 $n$ 列中选出 $m$ 个线性无关的列，构成 $m \times m$ 非奇异子矩阵 $B$。
- $B$ 的列向量 = 基向量
- $B$ 对应的变量 = 基变量（$m$ 个）
- 其余变量 = 非基变量（$n-m$ 个）
基解：令所有非基变量 = 0，解 $BX_B = b$，得到 $X_B = B^{-1}b$。再加上非基变量（都是 0），就构成一个基解。
基可行解：满足非负约束 $X_B \ge 0$ 的基解。
可行域顶点与基可行解的一一对应（定理 2）：
$X$ 是可行域顶点 $\iff$ $X$ 是基可行解。
最优解的存在性（定理 3）：
若线性规划问题有最优解，则一定存在一个基可行解是最优解。
这意味着：我们只需要在 有限个顶点（基可行解） 中寻找最优解！这正是单纯形法的理论基础。
几何概念与代数概念的对应：约束超平面的交点 → 基解；可行域的顶点 → 基可行解。
单纯形法：完整算法
单纯形法的核心思想：从一个基可行解出发，沿着可行域的边走向相邻的基可行解，每一步都使目标函数值改善，直到找到最优解。
三个关键步骤：
第一步：初始基可行解的确定
利用松弛变量构造单位阵作为初始可行基。因为每个 $\le$ 约束都加了一个松弛变量，这些松弛变量的系数列向量恰好构成一个 $m \times m$ 的单位阵 $I_m$！
令所有原始变量（非基变量）= 0，松弛变量 = $b_i$，就得到了初始基可行解。
对于 $\ge$ 和 $=$ 约束（没有现成的单位向量），后面会讲人工变量法（大M法）。
第二步：最优性判断——检验数
对任一基可行解 $X^{(0)}$，将目标函数用非基变量表示：
$$z = z^{(0)} + \sum_{j \in I_N} \sigma_j x_j$$
其中 $\sigma_j = c_j - \sum_{i=1}^{m} c_i a_{ij}$ 称为检验数（即单纯形表中的 $C_j - Z_j$）：
- 若所有 $\sigma_j \le 0$（max 问题），则当前解就是最优解。
- 若有 $\sigma_j > 0$，引入 $x_j$ 还能进一步改善目标值。
第三步：基可行解的迭代
入基变量：选择 $\sigma_k = \max{\sigma_j \mid \sigma_j > 0}$，$x_k$ 成为新的基变量。
出基变量：计算 $\theta = \min\limits_{i}{ \frac{x_i^{(0)}}{a_{ik}} ;|; a_{ik} > 0 }$，对应的第 $l$ 个基变量出基。
如果所有 $a_{ik} \le 0$，则 $\theta$ 无法确定，意味着 $x_k$ 可以无限增大——无界解！
为什么选最小的 $\theta$？ 因为 $x_k$ 从 0 开始增大时，$\theta$ 是第一个被”挤到 0”的基变量所对应的值。如果 $x_k$ 超过 $\theta$，某个基变量就会变成负数，破坏可行性。
旋转运算（Pivot）：以 $a_{lk}$ 为主元素，做行初等变换，使 $x_k$ 列变成单位向量（第 $l$ 行为 1，其余行为 0），同时更新 $b$ 列。
目标函数值的变化：$z^{(1)} = z^{(0)} + \theta \cdot \sigma_k$，因为 $\theta > 0, \sigma_k > 0$，所以每次迭代目标函数都会严格改善（除非遇到退化）。
原有图表演示
课堂上的三个表格展示的就是上述过程：
特殊情况与解的判别
综合 PPT 和运筹学理论，线性规划解的判别标准如下：
| 条件 | 结论 |
|---|---|
| 所有 $\sigma_j \le 0$，且非基变量的 $\sigma_j < 0$ | 唯一最优解 |
| 所有 $\sigma_j \le 0$，存在某个非基变量的 $\sigma_j = 0$ | 无穷多最优解 |
| 存在 $\sigma_j > 0$，且该列所有 $a_{ij} \le 0$ | 无界解（目标值可无限改善） |
| 所有 $\sigma_j \le 0$，但有人工变量在基中不为 0 | 无可行解 |
判断流程图：
1 | 计算所有 σj |
退化与循环
退化：某个基变量取值为 0 的现象。此时 $\theta = 0$，目标函数在这一步不会改善（$z^{(1)} = z^{(0)} + 0 \cdot \sigma_k = z^{(0)}$）。
退化可能导致循环：经过若干次迭代后又回到同一个基可行解。
对策：
- Bland 规则：入基选索引最小的正检验数对应的变量；出基在 $\theta$ 相同时选索引最小的基变量
- 摄动法：对 $b$ 列做微小扰动打破退化
- 在实际计算中，退化导致循环的情况比较罕见
人工变量法（大M法）
对于含 $\ge$ 约束或 $=$ 约束的问题，标准化后没有现成的单位矩阵作为初始基。
大M法的思想：在 $\ge$ 和 $=$ 约束中引入人工变量 $x_{ai} \ge 0$，同时在目标函数中加一个巨大的惩罚项 $-M \cdot x_{ai}$（max 问题，$M \to +\infty$）。如果最优解存在，人工变量一定会被逼到 0。
转化规则：
- $\sum a_{ij}x_j = b_i$ → $\sum a_{ij}x_j + x_{ai} = b_i$
- $\sum a_{ij}x_j \ge b_i$ → $\sum a_{ij}x_j - x_{si} + x_{ai} = b_i$（减剩余变量后，再加人工变量）
- 目标函数（max）：$z = \sum c_jx_j - M\sum x_{ai}$
例题：$\min z = -3x_1 + x_2$，s.t. $x_1 + x_2 \ge 2$, $x_1 - x_2 \le 1$, $x_1, x_2 \ge 0$
标准化：$\max z’ = 3x_1 - x_2$，约束变为 $x_1 + x_2 - s_1 = 2$, $x_1 - x_2 + s_2 = 1$。
第一个等式中没有现成的基变量，引入人工变量 $a_1$：$x_1 + x_2 - s_1 + a_1 = 2$。
目标函数变为 $\max z’ = 3x_1 - x_2 + 0s_1 + 0s_2 - Ma_1$。
初始基为 $(a_1, s_2)$，构造单纯形表求解。如果最终 $a_1 = 0$（出基），就得到了原问题的最优解；如果 $a_1$ 始终为正而所有检验数已经 ≤ 0，说明原问题无可行解。
两阶段法是大M法的改进版本，将求解分为两阶段：第一阶段先最小化人工变量之和（判断可行性），第二阶段再优化原目标函数。在实际计算中更稳定，但思想相同。
电力系统中的应用
回到我们最初的问题——经济调度。当发电机成本是线性的（或分段线性的），ED 问题可以写成标准线性规划：
$$\min \sum_{i=1}^{N} c_i P_i$$
$$s.t. \quad \sum_{i=1}^{N} P_i = P_D, \quad P_i^{\min} \le P_i \le P_i^{\max}$$
这就是一个标准的 LP 问题！用单纯形法可以直接求解。
更有意思的是——还记得经济调度中的增量成本 $\lambda$ 吗？在单纯形法的最终表中，松弛变量的检验数的相反数就是影子价格（对偶变量）。对于约束 $P_1 + P_2 + \dots = P_D$，对应的松弛变量的影子价格，就是系统 $\lambda$——多家 1MW 负荷需要多花的钱。
所以，经济调度和线性规划本质上是同一个框架：拉格朗日乘子法在连续可微条件下的特殊形式，单纯形法在线性条件下的通用形式。
更进一步，最优潮流（Optimal Power Flow, OPF）可以看作是带有电网潮流约束的扩展经济调度问题，既可以建模为线性规划（DC-OPF，直流潮流近似），也可以建模为非线性规划（AC-OPF，交流潮流精确模型）。
例题
例题一
所以X1=50，X2=0，Z=100
例题二
Optimal Power Flow（最优潮流）5.19
从 ED/LP 到 OPF
前面我们讨论的经济调度（ED）只关心有功功率平衡——所有发电机出力之和必须等于总负荷加损耗，但不关心这些功率怎么在电网中流动。
而最优潮流（OPF）更进一步：在满足经济调度的基础上，还要考虑电网的物理约束——包括节点电压必须在允许范围内、线路潮流不能越限、变压器分接头位置、无功补偿设备状态等。
从 LP 的角度看，当 ED 中加入了网络约束（line flow constraints），它就从纯经济调度变成了 OPF 的雏形。
OPF 的数学表述
OPF 是一个典型的非线性约束优化问题：
三类变量
| 类型 | 含义 | 例子 |
|---|---|---|
| 控制变量 $u$ | 可以主动调节的量 | 发电有功出力 $P_G$、发电机电压 $V_G$、变压器变比、移相器角度、投切电容器/电抗器状态、HVDC/FACTS 控制量、切负荷量 |
| 状态变量 $x$ | 随控制变量变化的系统响应 | 各节点电压幅值 $V_m$（除发电机节点）、各节点电压相角 $\theta_m$（除平衡节点） |
| 参数 $p$ | 已知且固定不变的系统特征 | 网络拓扑、线路参数 $(R, X, B)$、发电机成本函数、发电机限值、负荷、潮流和电压限值 |
经典目标函数
$$\min \sum_{i} F_i(P_{Gi}) \quad \text{最小化总发电成本}$$
其他可能的目标函数：
- 最小化系统损耗（Loss Minimization）
- 最小化控制变量的变化量（减少调节成本）
约束条件
等式约束 — 每个节点的有功和无功功率必须平衡（即潮流方程）：
$$P_{Gm} - P_{Dm} = P_m(\theta, V) \quad \forall m$$
$$Q_{Gm} - Q_{Dm} = Q_m(\theta, V) \quad \forall m$$
- 左边：节点 $m$ 的净注入功率（发电 - 负荷）
- 右边：由网络方程计算出的净注入功率（取决于全网电压和相角）
- 潮流方程的本质：找到一组 $(\theta, V)$ 使得计算注入 = 实际注入
不等式约束：
- 控制变量限值：$P_{Gi}^{\min} \le P_{Gi} \le P_{Gi}^{\max}$
- 线路潮流限值：$|S_{ij}| \le S_{ij}^{\max}$
- 电压限值：$V_m^{\min} \le V_m \le V_m^{\max}$
OPF 紧凑形式：
$$\min_u f(x, u)$$
$$s.t. \quad g(x, u) = 0 \quad \text{（等式约束：潮流方程）}$$
$$h(x, u) \le 0 \quad \text{（不等式约束）}$$
电网建模
瞬时量基础（稳态、恒频交流系统）
$$v(t) = V_{\max} \cos(\omega t + \theta_v)$$
$$i(t) = I_{\max} \cos(\omega t + \theta_I)$$
$$\omega = 2\pi f \quad \text{（$f = 60$ Hz 北美标准，单位为弧度/秒）}$$
复功率与视在功率
$$\vec{S} = P + jQ$$
$$S = \sqrt{P^2 + Q^2} \quad \text{（视在功率）}$$
支路潮流方程（Branch Flow Equation）
输电线路用 $\pi$ 型等值模型 表示，连接节点 $m$ 和 $n$：
- 串联导纳：$y_{mn} = g_{mn} + jb_{mn} = \frac{1}{z_{mn}} = \frac{1}{r_{mn} + jx_{mn}}$
- 其中 $g_{mn} = \frac{r_{mn}}{r_{mn}^2 + x_{mn}^2}$，$b_{mn} = -\frac{x_{mn}}{r_{mn}^2 + x_{mn}^2}$
- 并联导纳（线路充电电容）：两端各有一半 $\frac{jb_{ch, mn}}{2}$
节点功率平衡方程
对任意节点 $m$，净注入功率等于流向所有相邻节点的功率之和：
$$P_{Gm} - P_{Dm} = \sum_{n \in S_m} P_{mn}$$
$$Q_{Gm} - Q_{Dm} = \sum_{n \in S_m} Q_{mn}$$
其中 $S_m$ 为与节点 $m$ 直接相连的所有节点集合。
功率损耗
线路 $m$-$n$ 上的有功损耗：
$$P_{L, mn} = P_{mn} + P_{nm} = (V_m^2 + V_n^2 - 2V_m V_n \cos\theta_{mn}) \cdot g_{mn}$$
推导中利用了 $\cos\theta_{nm} = \cos\theta_{mn}$ 和 $\sin\theta_{nm} = -\sin\theta_{mn}$，使得含 $b_{mn}$ 的项对消。
AC 潮流方程（完整非线性形式）
以两节点系统为例，线路参数 $y_{12} = -5 + j5$（导纳矩阵元素），潮流方程如下：
有功平衡
$$P_{G1} - P_{D1} = V_1^2 \cdot G_{11} + V_1 V_2(G_{12}\cos\theta_{12} + B_{12}\sin\theta_{12})$$
$$P_{G2} - P_{D2} = V_2^2 \cdot G_{22} + V_2 V_1(G_{21}\cos\theta_{21} + B_{21}\sin\theta_{21})$$
无功平衡
$$Q_{G1} - Q_{D1} = -V_1^2 \cdot B_{11} + V_1 V_2(G_{12}\sin\theta_{12} - B_{12}\cos\theta_{12})$$
$$Q_{G2} - Q_{D2} = -V_2^2 \cdot B_{22} + V_2 V_1(G_{21}\sin\theta_{21} - B_{21}\cos\theta_{21})$$
变量：$V_1, V_2, \theta_1, \theta_2, P_1, P_2, Q_1, Q_2$
潮流方程的特点：非线性方程式，需给定初始解，用牛顿-拉夫森法迭代求解，雅可比矩阵在每一步更新。
牛顿-拉夫森法求解潮流
求解思路：每一步计算当前功率失配量 $\Delta P, \Delta Q$，通过雅可比矩阵 $J$ 得出修正量 $\Delta\theta, \Delta V$：
$$\begin{bmatrix} \Delta P \ \Delta Q \end{bmatrix} = J \begin{bmatrix} \Delta\theta \ \Delta V \end{bmatrix}$$
$$\begin{bmatrix} \Delta\theta \ \Delta V \end{bmatrix} = J^{-1} \begin{bmatrix} \Delta P \ \Delta Q \end{bmatrix}$$
然后更新 $\theta \leftarrow \theta + \Delta\theta$，$V \leftarrow V + \Delta V$，重复直到失配量足够小。
Jacobian 矩阵的典型元素（以 $\frac{\partial Q_1}{\partial V_1}$ 为例）：
$$Q_1 = -V_1^2 B_{11} + V_1 V_2(G_{12}\sin\theta_{12} - B_{12}\cos\theta_{12})$$
$$\frac{\partial Q_1}{\partial V_1} = -2V_1 B_{11} + V_2(G_{12}\sin\theta_{12} - B_{12}\cos\theta_{12})$$
$$\frac{\partial Q_1}{\partial V_2} = V_1(G_{12}\sin\theta_{12} - B_{12}\cos\theta_{12})$$
节点类型分类
| 节点类型 | 已知量 | 待求量 | 说明 |
|---|---|---|---|
| PQ 节点（负荷节点） | $P, Q$ | $V, \theta$ | 无无功功率控制能力 |
| PV 节点（发电机节点 / $V\theta$ 节点） | $P, V$ | $Q, \theta$ | 有无功功率控制能力 |
| Slack 节点（平衡节点） | $V, \theta$（$\theta=0$ 为参考角） | $P, Q$ | 补偿系统损耗，保证功率平衡 |
DC 潮流模型（线性化）
AC 潮流方程为非线性，难直接嵌入 LP。DC 潮流是一种线性近似：
关键假设：
- 线路电抗远大于电阻：$r_{mn} \ll x_{mn}$，故 $g_{mn} \approx 0$
- 节点电压幅值全为 1.0 pu
- 相角差很小：$\sin\theta_{mn} \approx \theta_{mn}$，$\cos\theta_{mn} \approx 1$
在此假设下，DC 潮流简化为：
$$P_{inj} = B \cdot \theta$$
其中 $B$ 是电纳矩阵（$B_{mm} = \sum_n \frac{1}{x_{mn}}$，$B_{mn} = -\frac{1}{x_{mn}}$）。
消除参考节点：设定一个参考节点（slack bus）的 $\theta_{\text{ref}} = 0$，从 $B$ 矩阵中删去其对应的行和列，得到降阶系统的 $\theta$。
转移因子（Shift Factor, SF）
SF 的含义：某个节点的净注入功率变化 1MW 时，会在各条线路上引起多少潮流变化。
$$SF = X^{-1} K_L^T BI$$
其中：
- $X$：支路电抗矩阵（对角阵）
- $K_L$：节点-支路关联矩阵（bus-branch incidence matrix）
- $BI$：$B$ 矩阵的”逆”（降阶 $B$ 矩阵的逆）
线路潮流 = SF × 节点净注入：
$$\begin{bmatrix} PL_{1 \to 2} \ PL_{1 \to 3} \ PL_{2 \to 3} \end{bmatrix} = [SF] \begin{bmatrix} P_{G1} - P_{D1} \ P_{G2} - P_{D2} \ P_{G3} - P_{D3} \end{bmatrix}$$
每条线路的潮流就是 SF 矩阵对应行与净注入向量的线性组合。
考虑网络约束的经济调度（SCED）
当我们将网络约束加入经济调度，就得到了安全约束经济调度（Security-Constrained Economic Dispatch, SCED）。
核心思想：先解 ED 得到发电计划，再跑一次潮流检查线路是否越限，只对越限的线路加入约束，避免一开始就把所有线路约束全加进去。
LP 形式的网络约束经济调度：
$$\min ; (1500 + 7PX_1) + (4000 + 8PX_2) + 10PX_3$$
$$s.t. \quad (100 + PX_1) + (200 + PX_2) + PX_3 = 600 \quad \text{（功率平衡）}$$
$$0 \le PX_1 \le 200, \quad 0 \le PX_2 \le 200, \quad 0 \le PX_3 \le 60$$
$$0.3PX_1 + 0.4PX_2 + 0.2PX_3 \le 90 \quad \text{（线路潮流约束）}$$
其中 $PX_i$ 是发电机在最小出力以上的增量出力，$P_i = P_i^{\min} + PX_i$。成本函数 $F_i = F_{i0} + c_i \cdot PX_i$，$F_{i0}$ 为最小出力处的成本（沉没成本）。
SF 矩阵中对应行的系数就是 LP 约束中 $PX_i$ 的系数。
损耗最小化 — 迭代 LP 方法
线路损耗最小化是一个非线性问题，但可以通过迭代求解 LP 逼近：
- 在当前运行点（由潮流计算得到）线性化损耗函数
- 求解 LP 得到 $\Delta V$（电压幅值修正）和 $\Delta Q$（无功出力修正）
- 更新 $V \leftarrow V + \Delta V$，$Q \leftarrow Q + \Delta Q$
- 重新求解潮流（非线性 AC 潮流），得到新的运行点
- 回到步骤 1，重复迭代直到收敛
损耗最小化是一个迭代过程（iterative process），每一步 LP 只是近似的线性化。
电压幅值的可调范围（约束）：
$$\Delta V_{m, \min} = V_{m, \min} - V_{m, 0}$$
$$\Delta V_{m, \max} = V_{m, \max} - V_{m, 0}$$
2 节点系统算例
一个 2 节点系统，$S_{g,1} = 0.7301 + j0.8301$，$S_{D1} = 0.6 + j0.3$，$S_{D2} = 0.1 + j0.5$（均 pu），线路 $Z = 0.1 + j0.1$ pu，已知 $V_1 = 0.9954\angle 0^\circ$，$V_2 = 0.93\angle -4.76^\circ$，系统损耗 $P_L = 0.03006$ pu。
约束：
- $0.9 \le V_1 \le 1.1$，$0.8 \le V_2 \le 0.93$
- $0.7 \le Q_{G1} \le 1.1$
用迭代 LP 进行损耗最小化求解：
LP 问题形式（MATLAB linprog 风格）：
- 决策变量：$x = [\Delta V_1, \Delta V_2, \Delta Q_1, \Delta Q_2]^T$
- 目标：$f = [0.751338, -0.72634, 0, 0]$
- 等式约束 $A_{eq} \cdot x = b_{eq}$，其中 $b_{eq} = [0, 0]^T$
- 变量上下界 $lb \le x \le ub$
一次 LP 求解得到 $\Delta V$ 后更新 $V_1 = V_1 + \Delta V_1 = 0.9843$，且 $fval = -0.0122$ 表示损耗减少了——$\Delta P_{Loss} = -0.0122$ pu。
State Estimation & FDI Attacks（状态估计与虚假数据攻击）5.19
从 OPF 到状态估计
OPF 的求解依赖于准确的系统参数和测量数据。然而实际中的测量值总是带有误差。状态估计（State Estimation）的目的就是：利用带有误差的冗余测量数据，估计出最可能的真实系统状态。
状态估计基本模型
测量模型
$$z = h(x) + e$$
- $z$：实际测量值向量（已知）
- $x$：真实系统状态（未知，需要估计）
- $h(x)$：测量函数（描述状态如何映射为测量量）
- $e$：测量误差
状态估计问题
残差 = 实际测量值 − 测量函数计算值：
$$r = z - h(x)$$
状态估计的目标就是：找到使某个目标函数最小化的 $\hat{x}$。
加权最小二乘（WLS）估计
假设测量误差 $e_i$ 相互独立、服从零均值正态分布，权重取误差方差的倒数：
$$\min J(x) = \sum_i \frac{(z_i - h_i(x))^2}{\sigma_i^2}$$
直流模型下的线性状态估计
对于 DC 潮流模型，$h(x) = Hx$（$H$ 是常数矩阵），状态估计有闭式解（closed-form solution）：
$$\hat{x} = (H^T R^{-1} H)^{-1} H^T R^{-1} z$$
其中 $R$ 为测量误差的协方差矩阵（对角阵，对角线元素为各测量误差的方差 $\sigma_i^2$）。
注意：这里假设所有仪表都有小误差（正常测量环境）。WLS 对异常值敏感，需要配以坏数据检测。
坏数据检测（Bad Data Detection）
残差与误差的关系
$$\text{残差 } r = S \cdot e \quad \text{（$S$ 是残差灵敏度矩阵，奇异矩阵）}$$
- 残差 $r$ 是误差 $e$ 的线性变换
- $S$ 是奇异的（non-invertible），意味着某些攻击组合可以”藏”在 $S$ 的零空间中
归一化残差
$$\text{归一化残差} = \frac{r_i}{\sqrt{\Omega_{ii}}}$$
其中 $\Omega$ 是残差的协方差矩阵。
最大归一化残差策略：
- 如果测量集中存在单个坏数据，最大归一化残差将对应该坏测量
- 若假设 $r$ 服从标准正态分布 $N(0,1)$，则当 $|r_i| > 3$ 时，我们有充分理由怀疑测量 $i$ 是坏数据
卡方检验（Chi-square Test）
卡方分布：$K$ 个独立标准正态随机变量的平方和的分布，自由度为 $K = m - n$（$m$ 为测量数，$n$ 为状态变量数）。
统计量 $J(\hat{x}) \sim \chi^2(m-n)$。
$$P{J(\hat{x}) > \chi^2_{\alpha}(m-n)} = \alpha$$
- $\chi^2_{\alpha}(m-n)$：临界值（由显著性水平 $\alpha$ 和自由度决定）
- $\alpha$：显著性水平
- $1 - \alpha$：置信度
若 $J(\hat{x})$ 超过临界值，则判定存在坏数据。
虚假数据注入攻击（False Data Injection Attacks, FDI）
核心问题
坏数据检测依赖于残差。如果攻击使得残差保持不变，则传统检测手段失效。
FDI 攻击通过协同篡改多个测量值，使得：
- 系统状态被错误估计（可被导向攻击者期望的任意值）
- 残差不变 → 传统坏数据检测无法检测
攻击条件：攻击者必须掌握 $H$ 矩阵（网络拓扑 + 线路参数）。
两种典型的 FDI 攻击
1. 负荷重分配攻击（Load Redistribution Attack, LRA）
- 原理：篡改不同节点的负荷测量值，使得总负荷量保持不变，但负荷分布发生了变化
- 假设：对于大型发电机，其出力测量不能被攻击；零注入节点的测量也不能被攻击
- 负荷测量可在一定范围内被篡改（变化百分比受限制）
- 篡改负荷测量后，相应的支路潮流测量也必须协同篡改，保证残差不变
双层优化问题（Bilevel Programming）：
- 上层问题（攻击者模型）：在限制范围内篡改负荷测量，最大化系统调度成本
- 下层问题（系统操作员的调度模型）：在收到虚假数据后，按最优调度（OPF/ED）运行
- 攻击者知道操作员的调度策略（H 矩阵已知），可以预测操作员的反应
2. 拓扑保持攻击（Topology Preserving Attack, TPA）
- 场景：某条线路因物理故障或网络事件而实际中断
- 攻击目标：伪造测量数据，使状态估计结果认为该线路仍在正常运行
- 攻击后果：操作员看不到线路中断，可能导致连锁故障
- 这比 LRA 更危险，因为物理拓扑已变而操作员毫不知情
FDI 攻击小结
| 攻击类型 | 篡改对象 | 攻击效果 | 约束 |
|---|---|---|---|
| LRA | 负荷测量 | 负荷看似被重新分配，总负荷不变；导致操作员采取次优调度 | 发电机出力和零注入节点不可攻破；负荷篡改范围有限 |
| TPA | 状态测量 | 掩盖线路中断，操作员认为拓扑正常 | 需要协同篡改大量相关测量值 |
关键结论：FDI 攻击能绕过传统坏数据检测的根本原因在于——攻击者利用了 $S$ 矩阵的零空间。只要攻击向量 $a$ 满足 $a = Hc$（$c$ 是任意非零向量），攻击就不会改变残差。
5.19 课程总结
本次课（5.19）涵盖了两大主题：
OPF 与 LP 在 OPF 中的应用：
- OPF 的完整数学建模：控制变量、状态变量、参数、目标函数、约束
- 电网建模：从瞬时量到复功率，从支路潮流到节点功率平衡
- AC 潮流全非线性方程与牛顿-拉夫森迭代求解
- DC 潮流的线性化推导 —— 将非线性潮流转为 $P = B\theta$
- SF 矩阵 —— 将节点注入映射到线路潮流
- 网络约束经济调度的 LP 形式
- 损耗最小化的迭代 LP 框架
状态估计与 FDI 攻击：
- 测量模型 $z = h(x) + e$ 与加权最小二乘估计
- 残差分析与坏数据检测（归一化残差、卡方检验）
- FDI 攻击（LRA、TPA）利用残差不变量绕开检测
- 双层优化模型：攻击者 vs 系统操作员
