企业官网建设流程全解析

1. 从“完美复现”说起：这个项目到底在解决什么问题

前一段时间，有个做配网调度的朋友问我，能不能用分布式的方式做机组间的经济分配，不要总靠调度主站的集中计算。他给的背景很典型：随着分布式光伏、储能、小水电这些“散户”资源大量接入，传统“中心节点统管全局”的调度模式越来越吃力——不是算不动，而是通信依赖太重、拓扑一旦变化就得重新整定，而且各机组的数据全要汇总到一处，对数据隐私和容错都不友好。

当时我正好在整理多智能体一致性算法相关的算例，就把“基于多智能体系统一致性算法的电力系统分布式经济调度策略”这套Matlab实现完整落了一遍。这个项目的核心价值，说白了就是用一套迭代公式，让每台机组只跟自己的“邻居”交换信息，最终所有机组的增量成本（也就是边际成本）收敛到同一个值，同时还能满足总负荷需求和机组出力上下限约束。这是个非常经典的自组织协同问题，也是分布式优化在电力系统里最典型的落地场景之一。

先说这个项目适合谁看：如果你正在学多智能体一致性算法，想找一个“能跑通、能改参、能出图”的Matlab算例；或者你在做微电网、主动配电网的分布式调度仿真，需要一套可复现的对比基准；又或者你只是想理解“增量成本一致性”到底是怎么从理论变成代码的，这篇文章都值得看完。我会从问题建模、算法设计、代码实现、仿真调参到常见坑位，完整拆一遍。

需要说明一点：正文中的案例与数据均为虚构设计，用于演示算法流程，所有机构和人物均不指向任何真实对象。

2. 经济调度问题的标准建模与求解难点

2.1 目标函数与约束：让总发电成本最小

先把问题形式化。假设系统里有 N 台可控机组，每台机组的发电成本通常用二次函数近似：

C_i(P_i) = a_i * P_i^2 + b_i * P_i + c_i

其中 P_i 是第 i 台机组的出力，a_i、b_i、c_i 是成本系数。二次函数这个假设不是拍脑袋定的，它来自火电机组热耗特性的实测拟合——煤耗随着出力增加加速上升，所以二次项系数 a_i 一定是正的。

分布式经济调度的目标函数是让所有机组的总成本最小：

min Σ C_i(P_i)

约束条件分两类：一是等式约束，全校的发电总量要等于总负荷需求 Σ P_i = Pd；二是不等式约束，每台机组出力要在上下限之内，P_min_i ≤ P_i ≤ P_max_i。

这个模型如果放在集中式框架里，拉格朗日乘子法或者二次规划求解器一把梭就完事了。但分布式场景下有个关键差异：没有人知道全局的总负荷 Pd 是多少，也没有人掌握所有机组的成本参数。每个节点只有自己的本地信息和邻居交换来的信息。所以问题变成了：设计一套迭代规则，让全局目标函数的最优解被“分布式地”找出来。

2.2 为什么用一致性算法而不是集中式优化

集中式优化最大的隐患在于单点失效。调度主站一旦出问题，整个系统就失去了决策中心。另外，随着机组数量增加，集中式的通信带宽和计算压力都会成为瓶颈。更现实的问题是，很多新能源电站的数据并不愿意全部上送，分布式方法天然保护了各单元的隐私。

一致性算法的思想非常朴素：每个智能体维护一个状态变量，不断用邻居状态来修正自己，最终让全体状态趋于一致。数学上可以证明，在连通的无向图拓扑下，只要迭代权重矩阵设计合理，状态值会收敛到初始值的平均值。这个“平均值”性质在分布式经济调度中非常值钱——它能让所有机组在不知道全局信息的情况下，对某个公共量达成共识。

用生活化的类比来理解：一群不认识路的人要在一个地点集合，没有人有对讲机，只能告诉身边几个人自己的位置，然后根据身边人的位置调整自己的方向。只要这个“互相传递位置”的网络是连通的，所有人最终都能聚到同一个点。一致性算法干的就是这件事。

3. 分布式经济调度的核心思路：增量成本一致性

3.1 最优性条件：从拉格朗日函数到等微增率准则

要让总成本最小且满足供需平衡，经典做法是构造拉格朗日函数：

L = Σ C_i(P_i) - λ(Σ P_i - Pd)

对 P_i 求偏导并令其为零，得到：

dC_i/dP_i = λ

也就是所有机组的增量成本相等。这个条件在电力系统里叫“等微增率准则”——每台机组每多发一度电所增加的成本相等时，总成本达到最优。λ 的含义是系统的边际成本，单位是元/MWh。

这个条件太重要了。它把复杂的全局优化问题转化成了一个公共变量 λ 的求解问题：只要让每台机组的增量成本都收敛到同一个 λ*，并且总的出力恰好等于负荷，那么每个机组的出力就是全局最优解。

问题是，λ* 到底是多少？集中式方法可以直接解 KKT 条件，但分布式方法不知道负荷总量，只能通过迭代去“逼近”这个最优值。

3.2 一致性变量选什么：增量成本比出力更合适

有些初学者会想：既然要调度出力，那不如让每台机组直接对 P_i 做一致性协商，最后收敛到一个相同的出力值。这显然是错的——不同机组成本特性不同，最优出力本来就不一样，强行让出力一致会导致成本严重偏离最优。

正确做法是让每台机组维护一个“增量成本”作为一致性变量，记为 λ_i。每轮迭代中，机组 i 按照以下规则更新：

λ_i(k+1) = λ_i(k) - ε * Σ_{j∈N_i} (λ_i(k) - λ_j(k)) - δ * (P_i(k) - Pd_i_local)

其中 N_i 是机组 i 的邻居集合，ε 是耦合增益，δ 是功率偏差反馈增益，Pd_i_local 是机组 i 本地估计的负荷需求。这里面的第二项是标准的一致性协商项，让相邻机组的增量成本互相靠拢；第三项是功率平衡校正项，让发电总量向负荷需求靠拢。实际工程中，第三项里的 Pd_i_local 也需要通过一个分布式平均一致性估计器来获得，因为没有一个节点知道全局负荷。

得到 λ_i 之后，如何映射回机组出力？对二次成本函数求导，可以直接反解：

P_i = (λ_i - b_i) / (2 * a_i)

这个 P_i 还要做限幅处理：如果低于 P_min_i 就取 P_min_i，高于 P_max_i 就取 P_max_i。限幅之后，真实的增量成本就不再等于 λ_i 了，所以算法对 λ_i 的修正会体现出约束的作用——这恰是分布式算法处理不等式约束的巧妙之处。

3.3 算法为什么能收敛：拉普拉斯矩阵与收敛速度

一致性迭代的收敛性分析绕不开图论。假设通信拓扑是一个无向连通图，邻接矩阵记为 A，度数矩阵记为 D，拉普拉斯矩阵 L = D - A。标准的连续一致性迭代可以写成矩阵形式：

λ(k+1) = λ(k) - ε * L * λ(k)

因为拉普拉斯矩阵有一个特征值为 0，对应的特征向量是全 1 向量，所以在 0 特征值方向上的分量不会衰减；而其他特征值都是正实数，对应方向上的分量以 (1 - ε*λ_i)^k 的速度衰减。只要 ε 满足 0 < ε < 2/λ_max(L)，系统就能收敛到全体 λ 的平均值。

这里有个重要直觉：第二小特征值 λ_2(L)，也就是代数连通度，决定了一步迭代能把不一致的分量压下去多少。拓扑越紧密，λ_2 越大，收敛越快。所以在做仿真时，如果发现收敛很慢，首先怀疑的不是迭代步长，而是通信拓扑的连通度。

4. Matlab代码实现：从矩阵构建到迭代循环

4.1 参数初始化与通信拓扑构建

写代码的第一步不是写迭代，而是把系统参数和通信拓扑先搭好。我这里用 5 台机组的系统作为演示，成本系数和出力上下限都通过结构体统一管理，方便后续改成任意规模。

%% 系统参数初始化 % 机组成本系数 C_i = a_i * P_i^2 + b_i * P_i + c_i a = [0.0025; 0.0030; 0.0020; 0.0035; 0.0028]; % 二次项系数 b = [2.00; 1.80; 2.20; 1.60; 2.10]; % 一次项系数 c = [30; 25; 35; 28; 32]; % 常数项系数 % 出力上下限 P_min = [50; 40; 60; 30; 45]; P_max = [300; 250; 350; 200; 280]; % 总负荷需求 Pd = 800;

通信拓扑用邻接矩阵定义。为了让结果有对比意义，我先给一个环形拓扑——每台机组只跟相邻两台通信，这种拓扑连通度低，收敛偏慢，适合展示算法的收敛过程。

% 通信拓扑：环形拓扑，5台机组首尾相接 adj = zeros(5, 5); adj(1, [2, 5]) = 1; adj(2, [1, 3]) = 1; adj(3, [2, 4]) = 1; adj(4, [3, 5]) = 1; adj(5, [4, 1]) = 1; % 无向图对称化 adj = adj + adj';

邻接矩阵定义完之后，度矩阵和拉普拉斯矩阵可以直接算出来：

deg = diag(sum(adj, 2)); L = deg - adj;

需要注意，这里度矩阵是每个节点的邻居数量，对无向图来说直接按列求和即可。如果拓扑是有向图，拉普拉斯矩阵的定义要换成有向图的版本，迭代的收敛性分析会更复杂，收敛条件不再只是连通性，还要求拓扑包含有向生成树。

4.2 一致性迭代主循环：增量成本更新与功率平衡

接下来是核心迭代循环。为了体现分布式特性，每个节点只能访问自己的成本参数、本地功率估计以及邻居发来的增量成本值。我在这里把“分布式”体现在代码结构上——虽然Matlab里写成了矩阵操作，但每个节点 i 的更新公式确实只依赖邻居信息。

%% 迭代参数 iter_max = 300; % 最大迭代次数 epsilon = 0.05; % 一致性耦合增益 delta = 0.10; % 功率平衡反馈增益 %% 初始化 lambda = zeros(5, 1); % 增量成本初始值 P = zeros(5, 1); % 出力初始值 P_est = Pd / 5 * ones(5, 1); % 本地负荷估计，初始均分 % 记录迭代过程 lambda_history = zeros(5, iter_max); P_history = zeros(5, iter_max); sumP_history = zeros(1, iter_max); for k = 1:iter_max % 步骤1: 根据当前增量成本计算出力 for i = 1:5 P(i) = (lambda(i) - b(i)) / (2 * a(i)); % 出力限幅 if P(i) < P_min(i) P(i) = P_min(i); elseif P(i) > P_max(i) P(i) = P_max(i); end end % 步骤2: 一致项更新 —— 每个节点只与邻居交换lambda lambda_new = lambda; for i = 1:5 neighbors = find(adj(i, :) == 1); consensus_sum = 0; for j = neighbors consensus_sum = consensus_sum + (lambda(i) - lambda(j)); end lambda_new(i) = lambda(i) - epsilon * consensus_sum * 2 - delta * (P(i) - P_est(i)); end lambda = lambda_new; % 步骤3: 本地负荷估计的平均一致性更新 P_est_new = P_est; for i = 1:5 neighbors = find(adj(i, :) == 1); sum_est = 0; for j = neighbors sum_est = sum_est + (P_est(i) - P_est(j)); end P_est_new(i) = P_est(i) - epsilon * sum_est * 2; end P_est = P_est_new; % 记录 lambda_history(:, k) = lambda; P_history(:, k) = P; sumP_history(k) = sum(P); end

在这个循环里，有两个细节值得展开。

第一，乘法因子 2 是怎么来的。因为无向图的邻接矩阵是对称的，度矩阵的值是邻居个数，拉普拉斯矩阵每一行对角元是度，非对角元是 -1。按照标准一致性迭代 λ_new = λ - εLλ 展开，第 i 行就是 λ_i - ε * ( deg_i * λ_i - Σ_{j∈N_i} λ_j )，也就是 λ_i - ε * Σ_{j∈N_i} (λ_i - λ_j)。这和我代码里写的其实是一致的。注意我没有直接用矩阵乘法 L*lambda，而是用邻接关系遍历累加，这样代码的可读性和“分布式感”更强，也更方便后续改成真正的分布式通信仿真。

第二，功率平衡项的符号方向。当本地出力大于本地负荷估计时，说明这台机组可能发多了，λ 应该往下调一点，所以用减号。反过来，如果出力小于负荷估计，λ 会上升，鼓励这台机组多发电。这个负反馈机制是算法能同时满足一致性条件和功率平衡的关键。

4.3 结果输出与误差分析

迭代结束后，先画三张图：第一张是所有机组的增量成本收敛过程，第二张是各机组出力曲线，第三张是总出力与负荷需求的偏差。我有一次跑完发现 λ 收敛到了同一个值，但总出力离 800 MW 差了 20 多 MW，问题就出在 δ 选得太小，功率平衡项修正力度不够。后来把 δ 调到 0.1，误差降到了 0.5 MW 以内。

还可以加一段计算收敛误差的代码，用来量化算法的收敛性能：

% 计算收敛误差: 每次迭代的lambda最大值与最小值之差 lambda_gap = zeros(1, iter_max); for k = 1:iter_max lambda_gap(k) = max(lambda_history(:, k)) - min(lambda_history(:, k)); end % 最终出力分配 P_final = P_history(:, end); % 总出力 P_total = sumP_history(end);

通过这个误差数组，可以观察收敛速度和稳态精度。理论上，在连续迭代下 λ 会严格收敛到同一个值，但由于我代码里加入了出力限幅的饱和环节，当某台机组触碰到出力边界时，未限幅的 λ 和限幅后的实际增量成本会不一致，这时迭代的收敛点不再是“所有 λ 相等”，而是“所有未饱和机组的 λ 相等，饱和机组的出力固定在边界上”。这是算法正确处理不等式约束时必然出现的行为，不算错误。

5. 仿真场景设计与参数调优实战

5.1 场景一：5机系统在环形拓扑下的收敛验证

先跑最基础的静止场景：负荷固定 800 MW，环形拓扑，观察增量成本和出力的收敛行为。我用 5 台机组的线性成本函数数据作为示例。算一下理论上的最优增量成本作为对照：不考虑限幅时，直接用集中式等微增率公式，总出力等于 800，可得 λ* ≈ 8.67。迭代结束后，五台机组的增量成本分别收敛到 8.67 附近，总出力也稳定在 800 MW。各机组出力分别为约 1334、1145、1618、1010、1273 MW（虚构示例，此处数值仅做形态示意，实际以运行结果为准），总成本比平均分配方案低了不少。

这里给个小复盘表格：

可以看到，一致性算法在几乎没有牺牲经济性的前提下，大幅降低了通信要求。

5.2 场景二：通信拓扑对收敛速度的影响

把环形拓扑换成全连接拓扑，也就是每台机组都能跟其他四台通信。因为拉普拉斯矩阵的第二小特征值变大了，收敛速度明显提升。我用同样的迭代步长 ε，环形拓扑大概需要 120 步左右把 λ 的差值压到 0.01 以内，全连接拓扑大约 40 步就够了。

这个观测对实际系统设计有直接指导意义：如果你的通信条件允许，适度增加每个节点的邻居数量，比单纯增大 λ 的迭代增益 ε 要更安全、更有效。因为增大 ε 容易诱发数值振荡甚至发散——我在后文的参数敏感性分析中会验证这一点。

5.3 场景三：负荷阶跃变化下的跟踪行为

在实际运行中，负荷不是恒定的。我设置在第 150 次迭代时，将总负荷从 800 MW 阶跃到 1000 MW，观察算法是否能跟踪到新的最优增量成本。结果是：负荷突变后，本地负荷估计 P_est 也会经过大约 30 次迭代收敛到新的平均值，λ 随之上升，出力总量逐步逼近 1000 MW。这个过程体现了算法对时变工况的适应性——它没有重新初始化，而是依靠持续的迭代和负反馈自行追踪。

需要提醒的是，负荷变化的速率不能太快，否则功率平衡项的修正跟不上。工程上常用的一种做法是给负荷估计器加一个低通滤波，让 P_est 的更新不要过于剧烈。

5.4 参数敏感性实验：ε 和 δ 怎么配

很多人拿到代码喜欢直接改 ε，发现改大了就发散，然后一头雾水。我做一个系统的参数扫描，在 0.01 到 0.2 范围内取不同的 ε 和 δ，记录 λ 到收敛所需要的迭代次数和最终功率偏差。

实验结果大致规律如下：ε 主导一致性收敛速度，但过大会导致振荡，因为迭代矩阵的谱半径超过了 1。δ 主导功率平衡的修正强度，太小会留下稳态偏差，太大会让 λ 的轨迹产生高频抖动。从经验来看，ε 建议取在 1/λ_max(L) 的一半左右，δ 取 ε 的 2 到 3 倍，二者按拓扑的规模做缩放。这种搭配在大多数小型测试系统里都能获得不错的收敛性能。

6. 完整代码框架与结果可视化

6.1 把分散的代码组织成可复跑的脚本

如果你只是跑一次看个结果，把上面的代码按顺序粘到一个 .m 文件里就够了。但如果你想把这个算例当成一个基础平台，后续要换拓扑、换成本参数、换负荷曲线，我建议花点时间整理成函数 + 脚本分离的结构。

一个清晰的文件组织可以是：

distributed_ed/ ├─ main_distributed_ed.m % 主脚本，包含参数设置、迭代循环、画图 ├─ build_topology.m % 函数：输入邻居关系，输出邻接矩阵、拉普拉斯矩阵 ├─ cost_params.m % 函数：返回成本系数和出力上下限 ├─ consensus_iteration.m % 函数：单步迭代逻辑，便于替换成不同算法变体 └─ plot_results.m % 函数：集中画图，便于统一风格

主脚本里只需要维护“算例配置”的代码块，迭代逻辑交给函数，这样你想对比不同拓扑对收敛速度的影响时，只需要改 build_topology 的输入，而不会碰坏核心算法代码。

6.2 核心函数：单步迭代的封装示例

function [lambda_new, P_est_new, P_out] = consensus_iteration(lambda, P_est, adj, a, b, P_min, P_max, epsilon, delta) N = length(lambda); lambda_new = lambda; P_est_new = P_est; % 根据当前 lambda 计算出力并限幅 P_out = zeros(N, 1); for i = 1:N P_out(i) = (lambda(i) - b(i)) / (2 * a(i)); P_out(i) = max(P_min(i), min(P_max(i), P_out(i))); end % 一致性协商 + 功率平衡修正 for i = 1:N neighbors = find(adj(i, :) == 1); for j = neighbors lambda_new(i) = lambda_new(i) - epsilon * (lambda(i) - lambda(j)); P_est_new(i) = P_est_new(i) - epsilon * (P_est(i) - P_est(j)); end lambda_new(i) = lambda_new(i) - delta * (P_out(i) - P_est(i)); end end

把迭代逻辑封装成函数的好处是，你可以很方便地替换成其它变体，比如加入投影算子的算法、带时延补偿的算法、或者自适应步长的算法。以后写论文做对比实验时，这种模块化结构能帮你节省大量时间。

6.3 结果可视化的几个关键视图

画图最忌讳的是信息堆在一起看不清。我习惯分四个子图展示结果：增量成本收敛曲线放在左上，各机组出力曲线放在右上，总出力与负荷对比放在左下，增量成本极值差（max-min）放在右下。前三个图让人直观感受收敛行为，第四个图给一个量化指标。

figure('Color', 'white', 'Position', [100, 100, 1000, 800]); subplot(2, 2, 1); plot(1:iter_max, lambda_history', 'LineWidth', 1.5); xlabel('迭代次数'); ylabel('增量成本 \lambda'); title('增量成本收敛过程'); legend('机组1', '机组2', '机组3', '机组4', '机组5', 'Location', 'best'); subplot(2, 2, 2); plot(1:iter_max, P_history', 'LineWidth', 1.5); xlabel('迭代次数'); ylabel('出力/MW'); title('各机组出力变化'); subplot(2, 2, 3); plot(1:iter_max, sumP_history, 'LineWidth', 2); yline(Pd, '--r', '负荷需求'); xlabel('迭代次数'); ylabel('总出力/MW'); title('总出力与负荷需求'); subplot(2, 2, 4); plot(1:iter_max, lambda_gap, 'LineWidth', 2); set(gca, 'YScale', 'log'); xlabel('迭代次数'); ylabel('增量成本差值'); title('一致性误差(对数坐标)');

用对数坐标看一致性误差是个很实用的习惯。线性坐标下误差降到 0.01 之后就“看起来”是零了，但换成对数坐标，能明显看到它有多个收敛阶段：先快速下降，然后进入平台，再缓慢下降。这种细节能帮你判断算法是不是卡在了某个环节，比如出力限幅导致的局部约束饱和。

7. 常见问题与排查技巧实录

7.1 迭代发散：ε 过大的典型表现

现象：λ 出现正负交替的锯齿形波动，数值越来越大，很快溢出。原因几乎都是 ε 超出了收敛边界 2/λ_max(L)。

排查步骤：先算出拉普拉斯矩阵的最大特征值，然后给 ε 一个远小于理论边界的值，比如理论边界的 1/4。不要直接猜，用特征值说话。你可以在代码里加一行调试输出：

eig_L = eig(L); epsilon_max = 2 / max(eig_L); fprintf('最大特征值: %.4f, 理论最大epsilon: %.4f\n', max(eig_L), epsilon_max);

有了这个数，你再设 ε 心里就有底了。另一个经验是：如果拓扑节点数增多，λ_max 也会增大，ε 的可取范围变窄，这时候与其继续调 ε，不如考虑改进拓扑结构，增连通度。

7.2 稳态功率偏差：δ 太小或负荷估计器收敛慢

现象：λ 已经相等了，但总出力和负荷之间差了十几兆瓦。这时候先看 P_est 是否已经收敛——把每轮迭代后的负荷估计误差打出来，如果 P_est 还没收敛，功率平衡修正项就在用错误的目标值做反馈。

解决办法有两个方向。一个是加大 δ，让功率修正更强；另一个是给 P_est 单独换一个更快的平均一致性增益，因为 P_est 的动态比 λ 慢的话，它会成为整个算法的瓶颈。我个人的调参顺序是：先保证 P_est 的收敛速度比 λ 快，再去调 ε 和 δ 的配合。

7.3 出力处于边界时增量成本不相等：不是 bug

这是最容易让人误判的情况。某台机组出力已经顶到上限，它的增量成本比别的机组低，但这是对的。因为如果它严格按照“增量成本相等”的原则，它的出力需要继续增加，但物理上不允许了。这时准确的说法是：这个机组处于“饱和”状态，它的 KKT 乘子参与平衡，但等式条件被不等式条件取代了。理解这一点对你判断算法是否正确至关重要。

我在第一次跑带限幅的仿真时就犯过这个迷糊，以为算法写错了，花了半天去查迭代公式。后来把饱和机组的实际增量成本曲线叠加到图上，才看明白：未饱和机组之间的 λ 是相等的，饱和机组的未受限 λ 则被“顶”在上面——这才是带约束最优的正确形态。

7.4 拓扑非连通：孤立节点无法收敛

如果你把通信拓扑改成了非连通的，比如两个孤立的子图，每个子图内部能达成“局部一致”，但整个系统的增量成本不会统一。这种问题在代码层面很隐蔽，因为矩阵运算不会报错，只是结果不对。建议每次换完拓扑后，顺手算一下拉普拉斯矩阵零特征值的重数——如果大于 1，说明图不连通。

zero_eig_count = sum(abs(eig(L)) < 1e-10); if zero_eig_count > 1 warning('通信拓扑不连通，存在孤立子图'); end

这个检查可以写成一个独立的脚本，在每次仿真前自动跑一遍，能省掉很多排查时间。

8. 实测经验与进一步扩展方向

整套代码跑通之后，我最大的感受是：分布式经济调度在原理上不难，难的是理解每个设计选择背后的权衡。比如为什么一致性变量选增量成本而非出力，为什么功率平衡项需要单独加一个负荷估计回路，为什么出力限幅会导致收敛点偏移。这些细节不通过实际跑代码很难真正内化。

我在实际调试中踩过最深的坑，是参数适配的问题。同一套 ε 和 δ，在 5 机系统里表现良好，换到 10 机系统之后又不行了。这说明参数必须跟着系统规模、拓扑结构和成本参数一起调。后来我做了一个简单的启发式缩放：ε 按 1/λ_max 的 0.2 倍取，δ 按 ε 的 2 倍取，基本上到任何规模的测试系统都能较快收敛。你可以把我的经验当成起点，再根据你自己的仿真场景微调。

后续这个平台还能往几个方向扩展。第一个方向是引入通信时延，模拟真实通信网络的不确定性。第二个方向是改成有向拓扑，研究非对称通信下的一致性收敛条件。第三个方向是把固定步长改成自适应步长，加快动态负荷下的响应速度。第四个方向是加入储能单元和可再生能源的随机出力，把模型扩展为多时段的动态经济调度。每个方向都有公开的经典文献可以参考，但在动手之前，先把眼前这套基本的一致性命题调到烂熟于心，比追新概念要重要得多。

如果要说一句最想告诉各位的实操心得，那就是：不要只盯着“收敛结果”看，多花时间把“不收敛”和“慢收敛”的原因弄清楚。把迭代曲线、特征值、拓扑结构三样对应起来看，收获比跑通一百组算例都大。