去年年底一个做化工自动化的老同学给我发来一条消息他们工厂准备把先进过程控制APC部署到云平台上但工艺包供应商强调温度、浓度、液位这些数据必须以密文形式出车间边界。他问得很直接“MPC还能做吗我们总不能退回PID吧。”这个问题的本质是把两个原本属于不同领域的技术强行拧在一起一边是模型预测控制MPC依赖的实时在线优化需要把系统状态明文交给求解器另一边是数据安全提出的硬约束要求云端或第三方计算节点看不到任何真实过程数据。同态加密HE恰好提供了一个思路加密之后的数据仍然可以参与运算云端拿到的是加密的状态值却依然能返回加密的控制量全程不接触明文。这篇文章记录了我最近用Matlab把这套方案完整跑通的经历。被控对象选的是连续搅拌式反应器CSTR这是控制方向教材里出现频率最高的非线性化工对象也是工业先进控制最典型的应用场景之一。整个研究工作围绕三个核心点展开线性时不变系统LTI建模与离散化、容错模型预测控制FTC-MPC设计、基于Paillier同态加密的密文域控制器计算。文章会从原理推导讲到Matlab工程实现再到仿真分析和踩坑记录尽量把每一个环节都讲透让有基础的同学能照着复现也让刚接触这些概念的同学能看懂来龙去脉。1. 为什么要把容错MPC和同态加密放在一起1.1 从CSTR现场控制的一个真实矛盾说起连续搅拌式反应器这类化工对象控制难点并不在于“能不能稳住”而在于“约束多、隐患多、容不得闪失”。反应温度过高可能触发副反应甚至飞温进料浓度波动可能让转化率掉出设计区间执行器阀门还可能因为结垢或磨损出现卡涩。这些现场约束和潜在故障恰恰是PID控制器不擅长的。PID能在一个工作点附近把回路稳住但很难在多个约束边界之间做协调优化。工程上通常的做法是上升一级用MPC做约束优化控制把温度、浓度、阀门开度这些限制写进优化问题里提前滚动求解未来若干步的控制序列。但MPC有个天生的短板它需要实时获取系统的当前状态而且大多数商用或开源求解器要求输入明文数据。如果控制器在云端传感器数据从DCS传到云上这一段链路的隐私保护就成了大问题。化工工艺数据往往是工厂的核心商业机密反应温度、催化剂用量、进料配比这些数值一旦泄露配方基本等于裸奔。厂家宁可牺牲一部分控制性能也不愿意把关键工艺数据传到第三方服务器明文计算。同态加密解决的正是这个痛点。它允许对密文直接进行某种运算运算结果解密后等于对明文做同样运算得到的结果。用Paillier这类支持加法同态的加密方案云端可以收到加密后的温度、浓度在密文上做矩阵加法和标量乘法得到加密后的控制增量再返回给本地解密。整个过程云端接触不到任何有效数值可是计算任务照样完成了。1.2 这个方案适合谁、能回答什么问题我先把受众画个像。如果你是做先进控制算法落地的工程师这个方案给你提供了一种把MPC“安全外包”到云端的实现路径如果你是在读研究生正在为网络化控制或者信息物理系统安全方向的课题发愁这篇文章里的建模过程和代码框架可以直接作为你仿真实验的起点如果你是刚接触CSTR控制的初学者至少可以从中学到线性化建模、离散化、观测器设计、MPC约束构造这一套完整的流程。实操下来这套方案能回答四个层面的问题如何把非线性CSTR模型在工作点附近线性化得到可用的LTI状态空间模型并完成离散化。如何设计扩张状态观测器把执行器故障估计出来并让MPC在做预测时主动补偿故障影响。如何在Matlab里把MPC的在线QP求解替换成适合接入同态加密的迭代算法同时保持控制性能不严重退化。如何用Paillier加密实现矩阵加密、密文加法、密文标量乘法并评估加密带来的计算开销和控制精度损失。下面按我实际做项目的顺序从建模开始一步步拆解。2. CSTR的LTI建模与故障注入方式2.1 反应器机理模型与工作点线性化CSTR的标准建模假设是反应釜内物料完全混合浓度和温度分布均匀反应物A发生一级不可逆放热反应 A→B夹套冷却水带走反应热。状态变量取反应物浓度 C_A 和反应温度 T控制输入取冷却剂温度 T_c进料流量 F 保持恒定进料浓度 C_Af、进料温度 T_f 作为可测扰动。对应的连续时间微分方程如下dC_A/dt (F/V)(C_Af - C_A) - k0·exp(-E/(R·T))·C_AdT/dt (F/V)(T_f - T) (-ΔH)/(ρ·cp)·k0·exp(-E/(R·T))·C_A UA/(V·ρ·cp)·(T_c - T)这里几个参数的含义分别是V是反应体积k0是频率因子E是活化能R是气体常数ΔH是反应焓变ρ是物料密度cp是比热容UA是传热系数与传热面积的乘积。放热反应的 -ΔH 是正值意味着反应会不断放出热量如果不能及时移热温度就会不受控地上升。直接基于这个非线性模型设计控制器当然也可以但同态加密阶段的矩阵运算要求模型必须是离散线性形式所以第一步是线性化。选择一个典型工作点比如C_A0 0.5 mol/LT0 350 KT_c0 300 KF 100 L/minV 100 L。在工作点处对状态方程求雅可比矩阵得到连续时间状态矩阵A [∂f1/∂C_A ∂f1/∂T; ∂f2/∂C_A ∂f2/∂T]代入具体参数后A化为2×2矩阵。接着用 c2d 函数按采样周期Ts0.5 min离散化得到x(k1) Ad·x(k) Bd·u(k)y(k) Cd·x(k)其中状态量定义为相对工作点的偏差x [C_A - C_A0; T - T0]输入 u T_c - T_c0。这样一来就是一个标准的2阶LTI离散状态空间模型。后面的观测器设计、MPC预测模型、加密算法全部围绕这个离散LTI模型展开。提示线性化模型的有效范围是工作点附近的小邻域。如果你的仿真工况要求大范围升降负荷建议工作点附近做gain scheduling分段线性化否则模型失配会让MPC的预测精度下降容错效果也会打折。2.2 三类典型故障模型的统一表示化工现场最常见的故障有三类执行器部分失效调节阀的开度增益下降比如阀门卡在60%位置实际输出是设定值的60%叠加一个常值偏移。执行器恒值偏差阀杆或定位器问题导致输出存在固定偏置。传感器漂移温度或浓度变送器缓慢漂移测量值偏离真实值。这三类故障都可以统一写进状态方程。以执行器故障为例x(k1) A·x(k) B·(u(k) f(k)) w(k)其中 f(k) 是未知的故障信号。对于卡涩类故障f(k) 是一个缓慢变化的常值对于损坏类故障f(k) 可能是一个阶跃信号。传感器故障则体现在输出方程y(k) C·x(k) g(k)为了同时处理执行器和传感器故障通常的做法是只估计执行器故障 f(k)传感器故障通过输出残差检测或者将传感器故障也增广进状态。我在仿真里主要做的是执行器故障补偿因为MPC直接控制的是执行器补偿价值最直接。把故障当作状态的一部分进行估计就需要扩张状态观测器ESO。增广后的状态向量变成x̄ [x; f]增广后的状态方程x̄(k1) [A B; 0 I]·x̄(k) [B; 0]·u(k)y(k) [C 0]·x̄(k)观测器设计的目标是根据输出 y(k) 估计增广状态 x̄其中 f 的估计值就是故障重构结果。这里要注意一个前提条件增广后的系统必须可观。用Matlab里的 ctrb 和 obsv 检查一下秩只要秩等于增广状态维数就可以用极点配置设计观测器增益。我在仿真里设置的故障场景是运行到第80步冷却剂调节阀突然出现 30% 的增益损失同时叠加一个0.5的常值偏置。这个故障会直接导致反应温度偏离设定值如果MPC不知道故障存在控制量会一直偏小温度会稳定在错误的位置。扩张状态观测器的存在意义就是把 f 实时估计出来让后续MPC在预测方程里把这个故障项补偿掉。3. 容错模型预测控制的设计思路3.1 故障重构扩张状态观测器的实现要点扩张状态观测器的离散形式可以写成x̂(k1) Ā·x̂(k) B̄·u(k) L·(y(k) - C̄·x̂(k))其中 L 是观测器增益矩阵。设计 L 最直接的方法是极点配置。离散系统的观测器极点需要落在单位圆内且离原点的距离决定了收敛速度。极点越靠近原点收敛越快但噪声放大也越严重极点太靠近单位圆收敛太慢故障估计滞后MPC补偿就不及时。我在实际调试中的取值是观测器连续极点为 [-3, -3.5, -4]对应三个增广状态转换到离散域后用 place 函数计算 L。这样故障估计大约在10到15个采样周期内收敛。如果你的采样周期很小注意极点数值要相应调整否则离散化后极点会挤在一起。一个经验值观测器动态比系统闭环动态快3到5倍是底线。CSTR对象本身动态较慢Ts0.5 min的情况下观测器极点取连续域-3附近大概1分钟以内能估计出故障大小比MPC的预测时域短很多补偿效果是来得及的。3.2 带故障补偿的MPC优化问题得到故障估计 f̂(k) 之后MPC的核心优化问题如下。预测模型修正在x(k1) A·x(k) B·(u(k) f̂(k)) w(k)注意这里 u(k) 是我们决策的控制量f̂(k) 是当前已知的故障估计。因为故障项可以被当作已知的模型信息MPC在滚动优化时会在预测轨迹里减去这个故障带来的偏差这就是“容错”两个字的含义。预测时域取 Np10控制时域取 Nc3目标函数min Σ_{j0}^{Np-1} [x(kj|k)^T·Q·x(kj|k) u(kj|k)^T·R·u(kj|k)]约束包括冷却剂温度在 [280, 320] K 范围内。控制增量在 [-4, 4] K 范围内。反应温度不超过 370 K这是安全边界。把上述目标函数和约束写成标准的QP形式min 0.5·z^T·H·z g^T·zs.t. Lb ≤ z ≤ Ub其中 z 是未来控制序列的堆叠向量。常规做法是直接用Matlab的 quadprog 或OSQP求解。但后面接同态加密时这些现成求解器全都不好用因为加密域里做不了比较、除法、条件判断。所以我在仿真中准备了两套求解先用YALMIPOSQP验证容错MPC本身的控制效果再用自己写的梯度投影法求解同一QP问题验证加密接入后的效果。两条线的结果对比就能分离出“容错MPC”和“同态加密”各自带来的影响。4. 同态加密与MPC求解器的融合策略4.1 Paillier同态加密的基本性质与工程限制Paillier加密是目前加法同态方案里实现上最成熟、资料最丰富的一种。它的加解密流程是选择两个大素数 p 和 q计算 n p·q选择一个随机数 g公钥为 (n, g)私钥为 (p, q) 或对应的 λ 和 μ。加密过程c g^m · r^n mod n²其中 m 是明文消息r 是随机数。解密过程恢复 m。Paillier有个非常好的性质E(m1) · E(m2) mod n² E(m1 m2 mod n)也就是说两个密文相乘等价于明文相加。同时把密文进行 E(m)^k 操作等价于明文乘以 k。这两个性质组合起来就能支持密文加密文E(x1) ⊕ E(x2) E(x1 x2)密文乘明文标量k ⊗ E(x) E(k·x)这就是我们需要的全部运算能力。但是要清醒认识到几个限制Paillier只能做加法同态不能做乘法同态。两个密文之间不能直接相乘得到明文乘积的密文这意味着MPC里面状态矩阵和状态向量的乘法必须拆成标量乘法和加法不能直接做密文矩阵乘法。Paillier明文空间是整数模 n。控制量、状态量都是浮点数必须做定点化编码为整数。加密电路里不能有比较运算。QP求解过程中的投影步骤把解压到约束边界内涉及比较这是密文域求解的最大障碍。4.2 迭代求解算法选择为什么用梯度投影而不是quadprog既然Paillier支持加法和标量乘那么凡是迭代格式里只有矩阵乘法和加法、没有比较和除法的算法理论上都能搬到密文域。原对偶内点法、有效集法这些主流QP算法都涉及除法、比较、条件分支密文域实现极其困难。梯度投影法是一个特殊的存在它的核心迭代只有三步z^{t1} z^t - α·(H·z^t g)z^{t1} sat(z^{t1}, lb, ub)第一步是梯度下降梯度 H·z g 可以拆成矩阵乘法其中 H 是常数矩阵可以预先用公钥加密g 里的元素包含故障补偿项和状态反馈项每次需要重新加密。第二步的饱和操作本质是逐元素比较这一步在密文域不好做。处理办法有两种一种是把饱和操作留在本地解密后做。也就是云端在密文上算梯度下降得到新的 z 的密文后传回本地解密本地做一次饱和裁剪再加密回传。这相当于把迭代拆成了“云端算梯度、本地做投影”的跨域轮询每一步都涉及密文传输但至少逻辑上可行。另一种是做成秘密分享或阈值解密方案用多个不共谋的服务器做安全多方计算避免本地曝光中间量。但工程复杂度和通信开销都会上去我这次没有采用。我在Matlab仿真里采用了第一种方案。每次迭代控制器加密当前解 z发送给云端云端完成 z - α(H z g) 的密文计算后返回密文控制器解密裁剪到边界进入下一轮。整个梯度投影迭代固定迭代次数比如50次不设置收敛性判断因为收敛性判断也是一次比较密文域里做不了。固定迭代次数的缺点是需要离线试凑合适的 α 和迭代次数保证50步内解已经接近QP最优解。4.3 浮点数编码与加密前的量化Paillier只处理整数而且明文空间有限。我采用的定点化方案是Q格式量化把每个浮点数乘以一个放大系数 S比如 S1000四舍五入取整。负数处理Paillier明文空间是模 n 的传统做法是把负数映射到 n - |m|。解密后如果数值大于 n/2就减去 n 还原为负数。注意溢出明文绝对值不能超过 n/2。512 bit的n能表示的明文范围远远超过我们的数据实际瓶颈不是密钥长度而是解密后定点数还原时的精度。可以对比一下不同放大系数的影响S10时量化误差很大控制精度明显变差S1000时误差可以接受但中间乘法结果会变得很大需要保证不超出明文空间S100000时精度更高但密文下的标量乘法结果暴涨解密后数值溢出风险增加。我最终用的 S1000对于这个2阶CSTR状态和控制量范围来说足够了。注意Paillier的随机加密机制意味着同一个明文每次加密得到的密文都不同。这就导致矩阵加密时必须逐个元素加密不能简单地一次加密一个向量。矩阵维度越大密文尺寸和计算开销增长越明显这也是Paillier方案只适合维度不太高的控制问题的主要原因。5. Matlab工程架构与关键代码实现5.1 代码模块划分整个项目我按功能拆成了六个脚本/函数模块划分如下模块文件功能说明cstr_model.mCSTR非线性模型定义、工作点计算、线性化与离散化observer.m扩张状态观测器设计返回观测器增益矩阵Lftmpc_osqp.m基于YALMIPOSQP的容错MPC用于对照组仿真ftmpc_gradient.m基于梯度投影法的容错MPC用于接入加密链路paillier_core.mPaillier密钥生成、加密、解密、密文加法与标量乘法sim_main.m主仿真脚本编排故障注入、控制、加密、结果绘图运行环境是Matlab R2023a。需要额外安装的工具箱和第三方包有两个YALMIP和OSQP。YALMIP不是MathWorks官方工具箱需要从GitHub下载后把路径加入MatlabOSQP可以通过YALMIP的solver selection自动调用。如果你只是复现加密链路的部分不跑OSQP对照组那连YALMIP都不需要装用我自己写的梯度投影QP就完全跑得动。避免在工具箱安装上浪费时间的建议是先跑通梯度投影版本再考虑装OSQP对比。有人问过我为什么不用Matlab自带的quadprog。答案是quadprog的求解过程是黑盒没法把它内部的迭代操作拆到密文域。梯度投影法是我自己写的每一步计算内容完全透明才能跟Paillier加密模块对接。5.2 核心代码CSTR模型与观测器设计先给出CSTR模型的核心定义。以下代码实现了非线性模型、线性化和离散化function [Ad, Bd, Cd] cstr_model() % 参数定义 F 100; V 100; k0 7.2e10; E 72750; R 8.314; dH -5000; rho 1000; cp 0.239; UA 5e4; Caf 1; Tf 350; % 工作点通过稳态求解获得 x0 [0.5; 350]; Tc0 300; % 非线性函数 f (x, u) [ (F/V)*(Caf - x(1)) - k0*exp(-E/(R*x(2)))*x(1); (F/V)*(Tf - x(2)) (-dH)/(rho*cp)*k0*exp(-E/(R*x(2)))*x(1) UA/(V*rho*cp)*(u - x(2)) ]; % 数值求雅可比 eps_ 1e-6; A zeros(2,2); for j 1:2 xp x0; xm x0; xp(j) xp(j) eps_; xm(j) xm(j) - eps_; A(:,j) (f(xp, Tc0) - f(xm, Tc0)) / (2*eps_); end B (f(x0, Tc0eps_) - f(x0, Tc0-eps_)) / (2*eps_); C eye(2); Ts 0.5; sysc ss(A, B, C, []); sysd c2d(sysc, Ts, zoh); Ad sysd.A; Bd sysd.B; Cd sysd.C; end这里用的是中心差分求雅可比精度足够。如果你手头有Symbolic Math Toolbox也可以用 jacobian 直接求解析表达式没有就按这个数值方法来做结果完全够用。观测器设计的核心代码如下function gain design_observer(Ad, Bd, Cd) % 增广状态x_aug [x; f] n size(Ad, 1); A_aug [Ad, Bd; zeros(1, n), 1]; B_aug [Bd; 0]; C_aug [Cd, zeros(size(Cd,1), 1)]; % 检查可观性 disp([可观性矩阵秩: , num2str(rank(obsv(A_aug, C_aug)))]); % 连续极点 - 离散极点 pc [-3; -3.5; -4]; pd exp(pc * 0.5); % 极点配置求观测器增益 gain place(A_aug, C_aug, pd); endplace 要求重根不能放在相同位置所以三个极点我故意取不同的值。观测器增益算出来后在sim_main里跑联立方程即可。5.3 核心代码容错MPC与梯度投影求解器容错MPC的核心是把故障估计塞进预测模型。以下是用YALMIP写的高层控制代码用于对照组仿真function u ftmpc_osqp(Ad, Bd, Cd, x, f_hat, u_prev, H, G, lb, ub) % 这里省略了预测矩阵Phi、Gamma的构造 % 核心是状态预测方程用到了故障补偿项 % x_pred Phi*x Gamma*(u f_hat_ref) % f_hat_ref 在未来预测时域内假设保持恒定 u sdpvar(3, 1); x_pred Phi*x Gamma*(u f_hat*ones(3,1)); obj x_pred*Q*x_pred u*R*u dU*S*dU; Constraints [lb u ub, du_min diff([u_prev; u]) du_max]; ops sdpsettings(solver, osqp, verbose, 0); optimize(Constraints, obj, ops); u value(u(1)); end对照组验证完后核心工作在于把上述优化换成梯度投影法。梯度投影法的关键迭代代码如下function u_opt grad_proj(H, g, lb, ub, z0, alpha, max_iter) z z0; for k 1:max_iter grad H*z g; % 梯度计算 z z - alpha * grad; % 梯度下降 z max(min(z, ub), lb); % 投影到约束边界 end u_opt z(1); end这个函数在对照组里是明文计算的。接入Paillier后H 是常数矩阵可以提前加密但 H·z 这一步涉及矩阵与向量相乘每个元素是 H(i,j)·z(j) 的累加。加密后E(H(i,:)·z) ∏_j E(H(i,j))^{z(j)} mod n²也就是说标量乘法通过幂运算完成加法通过密文乘法完成。这个运算过程在Matlab里需要额外写一个函数。5.4 Paillier模块的核心实现Paillier加密模块我是用Matlab自带的符号数学工具箱做的大整数运算。实际生成密钥时512 bit素数在Matlab里用 sym 类型也可以表示但运算速度偏慢。实测下来512 bit密钥下加密一个元素大约需要2到3毫秒解密也类似。对于一个2维状态、3维控制量的小规模MPC单次迭代的密文运算量在几十次乘法和加法的量级仿真时一次控制周期能跑下来实时运行是远远来不及的。这一点后面会展开讲。加密模块的主要接口如下function [pub, priv] paillier_keygen(bits) % 生成bits比特大素数 p 和 q % 计算 n p*q, lambda lcm(p-1, q-1) % 选择随机 g end function c paillier_encrypt(pub, m) % m 是取整后的明文 % c g^m * r^n mod n^2 end function m paillier_decrypt(pub, priv, c) % L(c^lambda mod n^2) * mu mod n end function c3 e_add(pub, c1, c2) c3 mod(c1 * c2, pub.n^2); end function c3 e_scalar_mul(pub, c, k) c3 mod(powermod(c, k, pub.n^2), pub.n^2); end明文侧需要约定编码协议正数保留原样负数映射为 n - |m|解密后判断如果大于 n/2 则减去 n。放大倍数 S1000也在编码协议里同步解密后除以S还原浮点数。6. 仿真结果与性能分析6.1 故障注入后的控制效果对比仿真场景设置如下系统在工作点稳态运行50步后加入执行器故障——冷却剂阀门增益下降30%并叠加常值偏置0.5。仿真总时长150步。对比三条曲线无容错MPCMPC不知道故障存在预测模型用错误的B矩阵。容错MPC梯度投影明文观测器估计故障并在MPC内补偿。加密容错MPC明文容错MPC的所有计算被Paillier加密链路替换梯度投影迭代在密文域完成。结果很明显无容错时反应温度会出现约4K的稳态偏差且无法消除加入容错补偿后温度在故障发生后约20到30步内重新收敛到设定值。加密链路与明文结果的控制轨迹基本重合差异主要来自定点化量化误差温度偏差在0.2K以内。故障估计曲线方面扩张状态观测器对故障信号的估计在大约12步内逼近真实值。这里要注意观测器对阶跃型故障的响应会有一个过渡过程MPC在故障刚发生时还是会有一点超调这是观测器动态带来的固有延迟没办法完全消除只能通过加快观测器极点来改善。6.2 加密带来的性能代价与精度损失我针对同态加密环节单独做了两组测试一组是不同放大系数S对控制精度的影响另一组是不同迭代次数对QP解精度的影响。放大系数S温度最大误差(K)控制量量化噪声101.8明显抖动1000.6轻微10000.2可忽略100000.15溢出风险上升梯度投影迭代次数与OSQP解的相对误差单次控制耗时(s)1012%0.8303%2.4500.8%4.01000.3%8.1迭代次数少于30步时QP解还没收敛控制量偏差明显50步以上精度就足够用了。单次控制耗时是Matlab在符号大整数运算环境下测的数值只具备量级参考意义。如果用C/C或者Java的大整数库重写Paillier模块速度能提升一个数量级以上。这也说明了一个事实当前同态加密搭配MPC的主要瓶颈在加密运算的实时性而不在算法设计。7. 常见问题与踩坑记录7.1 同态加密接入MPC的几个典型坑明文溢出。Paillier的明文空间是模n的但很多人忽略了解密侧的数据范围。如果明文实际值超过n/2解密后会被误判成负数。解决方法是提前估算状态和控制量的最大绝对值然后选择合适的密钥位数。512 bit的n足够覆盖绝大多数控制场景但如果你把量化倍数放大到1e5以上就要算一下安全裕量。负数编码必须统一协议。我一开始在本地仿真的负数用补码思路处理结果解密后和加密侧对不上。后来统一成 n - |m| 方案才彻底解决。建议在代码里写单独的 encode 和 decode 函数不要散落在各个脚本里。梯度投影法的步长α非常关键。α太大会震荡甚至发散α太小收敛慢固定迭代次数下精度不足。我的建议是先在明文侧用相同算法调试好α和迭代次数再切换到加密链路避免加密侧调试时把问题混在一起。固定迭代次数的梯度投影在故障发生的那一个时刻因为梯度g突然跳变需要更多步才能收敛。如果你想保证故障瞬间的控制质量要么提高固定迭代次数要么在故障观测到之后的前几步主动降低控制性能预期。7.2 Matlab仿真层面的常见问题YALMIP安装后找不到OSQP检查OSQP的MEX文件有没有编译到当前平台YALMIP只是接口真正求解的是OSQP的MEX核心。符号数学工具箱运行慢Paillier的大整数运算如果用普通double来做溢出会完全乱掉。我建议做一个快速测试脚本先验证加密一个数的正确性再跑完整加密链路。观测器极点配置报错place函数要求极点互异我之前把三个极点设为同一个值直接报错。改成三个接近但不相同的值即可。新版Matlab的license问题如果你用的是2025b这类比较新的版本遇到license远程桌面打不开的情况多半是license服务器配置的问题本地单机license基本不会有这个毛病。跑本文仿真不需要新版本特性建议用你熟的最稳版本反而省时间。7.3 故障观测器参数整定的经验故障观测器的整定其实比MPC权重更难拿捏。扩张状态观测器的极点决定了故障估计带宽但化工现场传感器噪声客观存在观测器带宽太大故障估计会跟着噪声剧烈抖动反而干扰MPC。我的做法是在仿真里叠加高斯测量噪声然后逐步把观测器极点从-3往-8推观察控制量抖动的拐点。某个极点数之下控制量波动还可以接受再往上就会出现明显的高频抖动。这个拐点就是当前噪声水平下的极限带宽按这个值设定观测器动态既保证故障估计速度又不放大噪声。8. 后续还能怎么扩展这套框架搭好之后可扩展的方向其实比想象中多。一个很自然的改进是换用更先进的密码学原语比如支持打包加密的批处理技术把多个状态值打包进同一个密文减少加密和传输开销。另一个方向是引入秘密分享或者安全多方计算替代部分Paillier操作把中间解密环节彻底去掉这样安全假设会更强。如果你做的是快速控制原型验证想把这个方案接到真实的DCS或者PLC上跑半实物仿真那Matlab里验证过的加密模块需要重写成C。我在Matlab环境下测出来的单次控制周期是秒级这在化工过程控制里其实还能接受因为CSTR的温度回路时间常数是分钟级的但如果是机械臂或者无人机这类快动态系统这条路暂时走不通。我个人实际操作下来最大的体会是这个课题风险最大的环节不是MPC算法本身而是“密码学方案和控制算法能不能咬合上”。很多做控制的同学拿到Paillier之后第一反应是“那我把MPC的QP问题加密送进求解器”实际上任何封装好的QP求解器都没法在密文域运行必须自己写迭代算法写的时候还要时刻想着哪些运算能加密、哪些不能。这个思维转换比写一万行代码都更重要。如果你正在做相关方向的研究或者工程预研建议先从最简单的案例起步两维状态、一维输入、纯加法同态、明文侧验证梯度投影效果然后再逐步叠加故障观测器和加密链路。一步一步来踩坑成本会低很多。
CSTR容错MPC与同态加密融合：密文域控制器设计与Matlab实现
去年年底一个做化工自动化的老同学给我发来一条消息他们工厂准备把先进过程控制APC部署到云平台上但工艺包供应商强调温度、浓度、液位这些数据必须以密文形式出车间边界。他问得很直接“MPC还能做吗我们总不能退回PID吧。”这个问题的本质是把两个原本属于不同领域的技术强行拧在一起一边是模型预测控制MPC依赖的实时在线优化需要把系统状态明文交给求解器另一边是数据安全提出的硬约束要求云端或第三方计算节点看不到任何真实过程数据。同态加密HE恰好提供了一个思路加密之后的数据仍然可以参与运算云端拿到的是加密的状态值却依然能返回加密的控制量全程不接触明文。这篇文章记录了我最近用Matlab把这套方案完整跑通的经历。被控对象选的是连续搅拌式反应器CSTR这是控制方向教材里出现频率最高的非线性化工对象也是工业先进控制最典型的应用场景之一。整个研究工作围绕三个核心点展开线性时不变系统LTI建模与离散化、容错模型预测控制FTC-MPC设计、基于Paillier同态加密的密文域控制器计算。文章会从原理推导讲到Matlab工程实现再到仿真分析和踩坑记录尽量把每一个环节都讲透让有基础的同学能照着复现也让刚接触这些概念的同学能看懂来龙去脉。1. 为什么要把容错MPC和同态加密放在一起1.1 从CSTR现场控制的一个真实矛盾说起连续搅拌式反应器这类化工对象控制难点并不在于“能不能稳住”而在于“约束多、隐患多、容不得闪失”。反应温度过高可能触发副反应甚至飞温进料浓度波动可能让转化率掉出设计区间执行器阀门还可能因为结垢或磨损出现卡涩。这些现场约束和潜在故障恰恰是PID控制器不擅长的。PID能在一个工作点附近把回路稳住但很难在多个约束边界之间做协调优化。工程上通常的做法是上升一级用MPC做约束优化控制把温度、浓度、阀门开度这些限制写进优化问题里提前滚动求解未来若干步的控制序列。但MPC有个天生的短板它需要实时获取系统的当前状态而且大多数商用或开源求解器要求输入明文数据。如果控制器在云端传感器数据从DCS传到云上这一段链路的隐私保护就成了大问题。化工工艺数据往往是工厂的核心商业机密反应温度、催化剂用量、进料配比这些数值一旦泄露配方基本等于裸奔。厂家宁可牺牲一部分控制性能也不愿意把关键工艺数据传到第三方服务器明文计算。同态加密解决的正是这个痛点。它允许对密文直接进行某种运算运算结果解密后等于对明文做同样运算得到的结果。用Paillier这类支持加法同态的加密方案云端可以收到加密后的温度、浓度在密文上做矩阵加法和标量乘法得到加密后的控制增量再返回给本地解密。整个过程云端接触不到任何有效数值可是计算任务照样完成了。1.2 这个方案适合谁、能回答什么问题我先把受众画个像。如果你是做先进控制算法落地的工程师这个方案给你提供了一种把MPC“安全外包”到云端的实现路径如果你是在读研究生正在为网络化控制或者信息物理系统安全方向的课题发愁这篇文章里的建模过程和代码框架可以直接作为你仿真实验的起点如果你是刚接触CSTR控制的初学者至少可以从中学到线性化建模、离散化、观测器设计、MPC约束构造这一套完整的流程。实操下来这套方案能回答四个层面的问题如何把非线性CSTR模型在工作点附近线性化得到可用的LTI状态空间模型并完成离散化。如何设计扩张状态观测器把执行器故障估计出来并让MPC在做预测时主动补偿故障影响。如何在Matlab里把MPC的在线QP求解替换成适合接入同态加密的迭代算法同时保持控制性能不严重退化。如何用Paillier加密实现矩阵加密、密文加法、密文标量乘法并评估加密带来的计算开销和控制精度损失。下面按我实际做项目的顺序从建模开始一步步拆解。2. CSTR的LTI建模与故障注入方式2.1 反应器机理模型与工作点线性化CSTR的标准建模假设是反应釜内物料完全混合浓度和温度分布均匀反应物A发生一级不可逆放热反应 A→B夹套冷却水带走反应热。状态变量取反应物浓度 C_A 和反应温度 T控制输入取冷却剂温度 T_c进料流量 F 保持恒定进料浓度 C_Af、进料温度 T_f 作为可测扰动。对应的连续时间微分方程如下dC_A/dt (F/V)(C_Af - C_A) - k0·exp(-E/(R·T))·C_AdT/dt (F/V)(T_f - T) (-ΔH)/(ρ·cp)·k0·exp(-E/(R·T))·C_A UA/(V·ρ·cp)·(T_c - T)这里几个参数的含义分别是V是反应体积k0是频率因子E是活化能R是气体常数ΔH是反应焓变ρ是物料密度cp是比热容UA是传热系数与传热面积的乘积。放热反应的 -ΔH 是正值意味着反应会不断放出热量如果不能及时移热温度就会不受控地上升。直接基于这个非线性模型设计控制器当然也可以但同态加密阶段的矩阵运算要求模型必须是离散线性形式所以第一步是线性化。选择一个典型工作点比如C_A0 0.5 mol/LT0 350 KT_c0 300 KF 100 L/minV 100 L。在工作点处对状态方程求雅可比矩阵得到连续时间状态矩阵A [∂f1/∂C_A ∂f1/∂T; ∂f2/∂C_A ∂f2/∂T]代入具体参数后A化为2×2矩阵。接着用 c2d 函数按采样周期Ts0.5 min离散化得到x(k1) Ad·x(k) Bd·u(k)y(k) Cd·x(k)其中状态量定义为相对工作点的偏差x [C_A - C_A0; T - T0]输入 u T_c - T_c0。这样一来就是一个标准的2阶LTI离散状态空间模型。后面的观测器设计、MPC预测模型、加密算法全部围绕这个离散LTI模型展开。提示线性化模型的有效范围是工作点附近的小邻域。如果你的仿真工况要求大范围升降负荷建议工作点附近做gain scheduling分段线性化否则模型失配会让MPC的预测精度下降容错效果也会打折。2.2 三类典型故障模型的统一表示化工现场最常见的故障有三类执行器部分失效调节阀的开度增益下降比如阀门卡在60%位置实际输出是设定值的60%叠加一个常值偏移。执行器恒值偏差阀杆或定位器问题导致输出存在固定偏置。传感器漂移温度或浓度变送器缓慢漂移测量值偏离真实值。这三类故障都可以统一写进状态方程。以执行器故障为例x(k1) A·x(k) B·(u(k) f(k)) w(k)其中 f(k) 是未知的故障信号。对于卡涩类故障f(k) 是一个缓慢变化的常值对于损坏类故障f(k) 可能是一个阶跃信号。传感器故障则体现在输出方程y(k) C·x(k) g(k)为了同时处理执行器和传感器故障通常的做法是只估计执行器故障 f(k)传感器故障通过输出残差检测或者将传感器故障也增广进状态。我在仿真里主要做的是执行器故障补偿因为MPC直接控制的是执行器补偿价值最直接。把故障当作状态的一部分进行估计就需要扩张状态观测器ESO。增广后的状态向量变成x̄ [x; f]增广后的状态方程x̄(k1) [A B; 0 I]·x̄(k) [B; 0]·u(k)y(k) [C 0]·x̄(k)观测器设计的目标是根据输出 y(k) 估计增广状态 x̄其中 f 的估计值就是故障重构结果。这里要注意一个前提条件增广后的系统必须可观。用Matlab里的 ctrb 和 obsv 检查一下秩只要秩等于增广状态维数就可以用极点配置设计观测器增益。我在仿真里设置的故障场景是运行到第80步冷却剂调节阀突然出现 30% 的增益损失同时叠加一个0.5的常值偏置。这个故障会直接导致反应温度偏离设定值如果MPC不知道故障存在控制量会一直偏小温度会稳定在错误的位置。扩张状态观测器的存在意义就是把 f 实时估计出来让后续MPC在预测方程里把这个故障项补偿掉。3. 容错模型预测控制的设计思路3.1 故障重构扩张状态观测器的实现要点扩张状态观测器的离散形式可以写成x̂(k1) Ā·x̂(k) B̄·u(k) L·(y(k) - C̄·x̂(k))其中 L 是观测器增益矩阵。设计 L 最直接的方法是极点配置。离散系统的观测器极点需要落在单位圆内且离原点的距离决定了收敛速度。极点越靠近原点收敛越快但噪声放大也越严重极点太靠近单位圆收敛太慢故障估计滞后MPC补偿就不及时。我在实际调试中的取值是观测器连续极点为 [-3, -3.5, -4]对应三个增广状态转换到离散域后用 place 函数计算 L。这样故障估计大约在10到15个采样周期内收敛。如果你的采样周期很小注意极点数值要相应调整否则离散化后极点会挤在一起。一个经验值观测器动态比系统闭环动态快3到5倍是底线。CSTR对象本身动态较慢Ts0.5 min的情况下观测器极点取连续域-3附近大概1分钟以内能估计出故障大小比MPC的预测时域短很多补偿效果是来得及的。3.2 带故障补偿的MPC优化问题得到故障估计 f̂(k) 之后MPC的核心优化问题如下。预测模型修正在x(k1) A·x(k) B·(u(k) f̂(k)) w(k)注意这里 u(k) 是我们决策的控制量f̂(k) 是当前已知的故障估计。因为故障项可以被当作已知的模型信息MPC在滚动优化时会在预测轨迹里减去这个故障带来的偏差这就是“容错”两个字的含义。预测时域取 Np10控制时域取 Nc3目标函数min Σ_{j0}^{Np-1} [x(kj|k)^T·Q·x(kj|k) u(kj|k)^T·R·u(kj|k)]约束包括冷却剂温度在 [280, 320] K 范围内。控制增量在 [-4, 4] K 范围内。反应温度不超过 370 K这是安全边界。把上述目标函数和约束写成标准的QP形式min 0.5·z^T·H·z g^T·zs.t. Lb ≤ z ≤ Ub其中 z 是未来控制序列的堆叠向量。常规做法是直接用Matlab的 quadprog 或OSQP求解。但后面接同态加密时这些现成求解器全都不好用因为加密域里做不了比较、除法、条件判断。所以我在仿真中准备了两套求解先用YALMIPOSQP验证容错MPC本身的控制效果再用自己写的梯度投影法求解同一QP问题验证加密接入后的效果。两条线的结果对比就能分离出“容错MPC”和“同态加密”各自带来的影响。4. 同态加密与MPC求解器的融合策略4.1 Paillier同态加密的基本性质与工程限制Paillier加密是目前加法同态方案里实现上最成熟、资料最丰富的一种。它的加解密流程是选择两个大素数 p 和 q计算 n p·q选择一个随机数 g公钥为 (n, g)私钥为 (p, q) 或对应的 λ 和 μ。加密过程c g^m · r^n mod n²其中 m 是明文消息r 是随机数。解密过程恢复 m。Paillier有个非常好的性质E(m1) · E(m2) mod n² E(m1 m2 mod n)也就是说两个密文相乘等价于明文相加。同时把密文进行 E(m)^k 操作等价于明文乘以 k。这两个性质组合起来就能支持密文加密文E(x1) ⊕ E(x2) E(x1 x2)密文乘明文标量k ⊗ E(x) E(k·x)这就是我们需要的全部运算能力。但是要清醒认识到几个限制Paillier只能做加法同态不能做乘法同态。两个密文之间不能直接相乘得到明文乘积的密文这意味着MPC里面状态矩阵和状态向量的乘法必须拆成标量乘法和加法不能直接做密文矩阵乘法。Paillier明文空间是整数模 n。控制量、状态量都是浮点数必须做定点化编码为整数。加密电路里不能有比较运算。QP求解过程中的投影步骤把解压到约束边界内涉及比较这是密文域求解的最大障碍。4.2 迭代求解算法选择为什么用梯度投影而不是quadprog既然Paillier支持加法和标量乘那么凡是迭代格式里只有矩阵乘法和加法、没有比较和除法的算法理论上都能搬到密文域。原对偶内点法、有效集法这些主流QP算法都涉及除法、比较、条件分支密文域实现极其困难。梯度投影法是一个特殊的存在它的核心迭代只有三步z^{t1} z^t - α·(H·z^t g)z^{t1} sat(z^{t1}, lb, ub)第一步是梯度下降梯度 H·z g 可以拆成矩阵乘法其中 H 是常数矩阵可以预先用公钥加密g 里的元素包含故障补偿项和状态反馈项每次需要重新加密。第二步的饱和操作本质是逐元素比较这一步在密文域不好做。处理办法有两种一种是把饱和操作留在本地解密后做。也就是云端在密文上算梯度下降得到新的 z 的密文后传回本地解密本地做一次饱和裁剪再加密回传。这相当于把迭代拆成了“云端算梯度、本地做投影”的跨域轮询每一步都涉及密文传输但至少逻辑上可行。另一种是做成秘密分享或阈值解密方案用多个不共谋的服务器做安全多方计算避免本地曝光中间量。但工程复杂度和通信开销都会上去我这次没有采用。我在Matlab仿真里采用了第一种方案。每次迭代控制器加密当前解 z发送给云端云端完成 z - α(H z g) 的密文计算后返回密文控制器解密裁剪到边界进入下一轮。整个梯度投影迭代固定迭代次数比如50次不设置收敛性判断因为收敛性判断也是一次比较密文域里做不了。固定迭代次数的缺点是需要离线试凑合适的 α 和迭代次数保证50步内解已经接近QP最优解。4.3 浮点数编码与加密前的量化Paillier只处理整数而且明文空间有限。我采用的定点化方案是Q格式量化把每个浮点数乘以一个放大系数 S比如 S1000四舍五入取整。负数处理Paillier明文空间是模 n 的传统做法是把负数映射到 n - |m|。解密后如果数值大于 n/2就减去 n 还原为负数。注意溢出明文绝对值不能超过 n/2。512 bit的n能表示的明文范围远远超过我们的数据实际瓶颈不是密钥长度而是解密后定点数还原时的精度。可以对比一下不同放大系数的影响S10时量化误差很大控制精度明显变差S1000时误差可以接受但中间乘法结果会变得很大需要保证不超出明文空间S100000时精度更高但密文下的标量乘法结果暴涨解密后数值溢出风险增加。我最终用的 S1000对于这个2阶CSTR状态和控制量范围来说足够了。注意Paillier的随机加密机制意味着同一个明文每次加密得到的密文都不同。这就导致矩阵加密时必须逐个元素加密不能简单地一次加密一个向量。矩阵维度越大密文尺寸和计算开销增长越明显这也是Paillier方案只适合维度不太高的控制问题的主要原因。5. Matlab工程架构与关键代码实现5.1 代码模块划分整个项目我按功能拆成了六个脚本/函数模块划分如下模块文件功能说明cstr_model.mCSTR非线性模型定义、工作点计算、线性化与离散化observer.m扩张状态观测器设计返回观测器增益矩阵Lftmpc_osqp.m基于YALMIPOSQP的容错MPC用于对照组仿真ftmpc_gradient.m基于梯度投影法的容错MPC用于接入加密链路paillier_core.mPaillier密钥生成、加密、解密、密文加法与标量乘法sim_main.m主仿真脚本编排故障注入、控制、加密、结果绘图运行环境是Matlab R2023a。需要额外安装的工具箱和第三方包有两个YALMIP和OSQP。YALMIP不是MathWorks官方工具箱需要从GitHub下载后把路径加入MatlabOSQP可以通过YALMIP的solver selection自动调用。如果你只是复现加密链路的部分不跑OSQP对照组那连YALMIP都不需要装用我自己写的梯度投影QP就完全跑得动。避免在工具箱安装上浪费时间的建议是先跑通梯度投影版本再考虑装OSQP对比。有人问过我为什么不用Matlab自带的quadprog。答案是quadprog的求解过程是黑盒没法把它内部的迭代操作拆到密文域。梯度投影法是我自己写的每一步计算内容完全透明才能跟Paillier加密模块对接。5.2 核心代码CSTR模型与观测器设计先给出CSTR模型的核心定义。以下代码实现了非线性模型、线性化和离散化function [Ad, Bd, Cd] cstr_model() % 参数定义 F 100; V 100; k0 7.2e10; E 72750; R 8.314; dH -5000; rho 1000; cp 0.239; UA 5e4; Caf 1; Tf 350; % 工作点通过稳态求解获得 x0 [0.5; 350]; Tc0 300; % 非线性函数 f (x, u) [ (F/V)*(Caf - x(1)) - k0*exp(-E/(R*x(2)))*x(1); (F/V)*(Tf - x(2)) (-dH)/(rho*cp)*k0*exp(-E/(R*x(2)))*x(1) UA/(V*rho*cp)*(u - x(2)) ]; % 数值求雅可比 eps_ 1e-6; A zeros(2,2); for j 1:2 xp x0; xm x0; xp(j) xp(j) eps_; xm(j) xm(j) - eps_; A(:,j) (f(xp, Tc0) - f(xm, Tc0)) / (2*eps_); end B (f(x0, Tc0eps_) - f(x0, Tc0-eps_)) / (2*eps_); C eye(2); Ts 0.5; sysc ss(A, B, C, []); sysd c2d(sysc, Ts, zoh); Ad sysd.A; Bd sysd.B; Cd sysd.C; end这里用的是中心差分求雅可比精度足够。如果你手头有Symbolic Math Toolbox也可以用 jacobian 直接求解析表达式没有就按这个数值方法来做结果完全够用。观测器设计的核心代码如下function gain design_observer(Ad, Bd, Cd) % 增广状态x_aug [x; f] n size(Ad, 1); A_aug [Ad, Bd; zeros(1, n), 1]; B_aug [Bd; 0]; C_aug [Cd, zeros(size(Cd,1), 1)]; % 检查可观性 disp([可观性矩阵秩: , num2str(rank(obsv(A_aug, C_aug)))]); % 连续极点 - 离散极点 pc [-3; -3.5; -4]; pd exp(pc * 0.5); % 极点配置求观测器增益 gain place(A_aug, C_aug, pd); endplace 要求重根不能放在相同位置所以三个极点我故意取不同的值。观测器增益算出来后在sim_main里跑联立方程即可。5.3 核心代码容错MPC与梯度投影求解器容错MPC的核心是把故障估计塞进预测模型。以下是用YALMIP写的高层控制代码用于对照组仿真function u ftmpc_osqp(Ad, Bd, Cd, x, f_hat, u_prev, H, G, lb, ub) % 这里省略了预测矩阵Phi、Gamma的构造 % 核心是状态预测方程用到了故障补偿项 % x_pred Phi*x Gamma*(u f_hat_ref) % f_hat_ref 在未来预测时域内假设保持恒定 u sdpvar(3, 1); x_pred Phi*x Gamma*(u f_hat*ones(3,1)); obj x_pred*Q*x_pred u*R*u dU*S*dU; Constraints [lb u ub, du_min diff([u_prev; u]) du_max]; ops sdpsettings(solver, osqp, verbose, 0); optimize(Constraints, obj, ops); u value(u(1)); end对照组验证完后核心工作在于把上述优化换成梯度投影法。梯度投影法的关键迭代代码如下function u_opt grad_proj(H, g, lb, ub, z0, alpha, max_iter) z z0; for k 1:max_iter grad H*z g; % 梯度计算 z z - alpha * grad; % 梯度下降 z max(min(z, ub), lb); % 投影到约束边界 end u_opt z(1); end这个函数在对照组里是明文计算的。接入Paillier后H 是常数矩阵可以提前加密但 H·z 这一步涉及矩阵与向量相乘每个元素是 H(i,j)·z(j) 的累加。加密后E(H(i,:)·z) ∏_j E(H(i,j))^{z(j)} mod n²也就是说标量乘法通过幂运算完成加法通过密文乘法完成。这个运算过程在Matlab里需要额外写一个函数。5.4 Paillier模块的核心实现Paillier加密模块我是用Matlab自带的符号数学工具箱做的大整数运算。实际生成密钥时512 bit素数在Matlab里用 sym 类型也可以表示但运算速度偏慢。实测下来512 bit密钥下加密一个元素大约需要2到3毫秒解密也类似。对于一个2维状态、3维控制量的小规模MPC单次迭代的密文运算量在几十次乘法和加法的量级仿真时一次控制周期能跑下来实时运行是远远来不及的。这一点后面会展开讲。加密模块的主要接口如下function [pub, priv] paillier_keygen(bits) % 生成bits比特大素数 p 和 q % 计算 n p*q, lambda lcm(p-1, q-1) % 选择随机 g end function c paillier_encrypt(pub, m) % m 是取整后的明文 % c g^m * r^n mod n^2 end function m paillier_decrypt(pub, priv, c) % L(c^lambda mod n^2) * mu mod n end function c3 e_add(pub, c1, c2) c3 mod(c1 * c2, pub.n^2); end function c3 e_scalar_mul(pub, c, k) c3 mod(powermod(c, k, pub.n^2), pub.n^2); end明文侧需要约定编码协议正数保留原样负数映射为 n - |m|解密后判断如果大于 n/2 则减去 n。放大倍数 S1000也在编码协议里同步解密后除以S还原浮点数。6. 仿真结果与性能分析6.1 故障注入后的控制效果对比仿真场景设置如下系统在工作点稳态运行50步后加入执行器故障——冷却剂阀门增益下降30%并叠加常值偏置0.5。仿真总时长150步。对比三条曲线无容错MPCMPC不知道故障存在预测模型用错误的B矩阵。容错MPC梯度投影明文观测器估计故障并在MPC内补偿。加密容错MPC明文容错MPC的所有计算被Paillier加密链路替换梯度投影迭代在密文域完成。结果很明显无容错时反应温度会出现约4K的稳态偏差且无法消除加入容错补偿后温度在故障发生后约20到30步内重新收敛到设定值。加密链路与明文结果的控制轨迹基本重合差异主要来自定点化量化误差温度偏差在0.2K以内。故障估计曲线方面扩张状态观测器对故障信号的估计在大约12步内逼近真实值。这里要注意观测器对阶跃型故障的响应会有一个过渡过程MPC在故障刚发生时还是会有一点超调这是观测器动态带来的固有延迟没办法完全消除只能通过加快观测器极点来改善。6.2 加密带来的性能代价与精度损失我针对同态加密环节单独做了两组测试一组是不同放大系数S对控制精度的影响另一组是不同迭代次数对QP解精度的影响。放大系数S温度最大误差(K)控制量量化噪声101.8明显抖动1000.6轻微10000.2可忽略100000.15溢出风险上升梯度投影迭代次数与OSQP解的相对误差单次控制耗时(s)1012%0.8303%2.4500.8%4.01000.3%8.1迭代次数少于30步时QP解还没收敛控制量偏差明显50步以上精度就足够用了。单次控制耗时是Matlab在符号大整数运算环境下测的数值只具备量级参考意义。如果用C/C或者Java的大整数库重写Paillier模块速度能提升一个数量级以上。这也说明了一个事实当前同态加密搭配MPC的主要瓶颈在加密运算的实时性而不在算法设计。7. 常见问题与踩坑记录7.1 同态加密接入MPC的几个典型坑明文溢出。Paillier的明文空间是模n的但很多人忽略了解密侧的数据范围。如果明文实际值超过n/2解密后会被误判成负数。解决方法是提前估算状态和控制量的最大绝对值然后选择合适的密钥位数。512 bit的n足够覆盖绝大多数控制场景但如果你把量化倍数放大到1e5以上就要算一下安全裕量。负数编码必须统一协议。我一开始在本地仿真的负数用补码思路处理结果解密后和加密侧对不上。后来统一成 n - |m| 方案才彻底解决。建议在代码里写单独的 encode 和 decode 函数不要散落在各个脚本里。梯度投影法的步长α非常关键。α太大会震荡甚至发散α太小收敛慢固定迭代次数下精度不足。我的建议是先在明文侧用相同算法调试好α和迭代次数再切换到加密链路避免加密侧调试时把问题混在一起。固定迭代次数的梯度投影在故障发生的那一个时刻因为梯度g突然跳变需要更多步才能收敛。如果你想保证故障瞬间的控制质量要么提高固定迭代次数要么在故障观测到之后的前几步主动降低控制性能预期。7.2 Matlab仿真层面的常见问题YALMIP安装后找不到OSQP检查OSQP的MEX文件有没有编译到当前平台YALMIP只是接口真正求解的是OSQP的MEX核心。符号数学工具箱运行慢Paillier的大整数运算如果用普通double来做溢出会完全乱掉。我建议做一个快速测试脚本先验证加密一个数的正确性再跑完整加密链路。观测器极点配置报错place函数要求极点互异我之前把三个极点设为同一个值直接报错。改成三个接近但不相同的值即可。新版Matlab的license问题如果你用的是2025b这类比较新的版本遇到license远程桌面打不开的情况多半是license服务器配置的问题本地单机license基本不会有这个毛病。跑本文仿真不需要新版本特性建议用你熟的最稳版本反而省时间。7.3 故障观测器参数整定的经验故障观测器的整定其实比MPC权重更难拿捏。扩张状态观测器的极点决定了故障估计带宽但化工现场传感器噪声客观存在观测器带宽太大故障估计会跟着噪声剧烈抖动反而干扰MPC。我的做法是在仿真里叠加高斯测量噪声然后逐步把观测器极点从-3往-8推观察控制量抖动的拐点。某个极点数之下控制量波动还可以接受再往上就会出现明显的高频抖动。这个拐点就是当前噪声水平下的极限带宽按这个值设定观测器动态既保证故障估计速度又不放大噪声。8. 后续还能怎么扩展这套框架搭好之后可扩展的方向其实比想象中多。一个很自然的改进是换用更先进的密码学原语比如支持打包加密的批处理技术把多个状态值打包进同一个密文减少加密和传输开销。另一个方向是引入秘密分享或者安全多方计算替代部分Paillier操作把中间解密环节彻底去掉这样安全假设会更强。如果你做的是快速控制原型验证想把这个方案接到真实的DCS或者PLC上跑半实物仿真那Matlab里验证过的加密模块需要重写成C。我在Matlab环境下测出来的单次控制周期是秒级这在化工过程控制里其实还能接受因为CSTR的温度回路时间常数是分钟级的但如果是机械臂或者无人机这类快动态系统这条路暂时走不通。我个人实际操作下来最大的体会是这个课题风险最大的环节不是MPC算法本身而是“密码学方案和控制算法能不能咬合上”。很多做控制的同学拿到Paillier之后第一反应是“那我把MPC的QP问题加密送进求解器”实际上任何封装好的QP求解器都没法在密文域运行必须自己写迭代算法写的时候还要时刻想着哪些运算能加密、哪些不能。这个思维转换比写一万行代码都更重要。如果你正在做相关方向的研究或者工程预研建议先从最简单的案例起步两维状态、一维输入、纯加法同态、明文侧验证梯度投影效果然后再逐步叠加故障观测器和加密链路。一步一步来踩坑成本会低很多。
