把控制器交给一个你管不到的计算环境还要让它继续做容错控制听起来就像让陌生人替你保管家门钥匙同时保证他永远不会偷看屋里放了什么。这正是我最近一直在折腾的问题线性时不变系统LTI的模型预测控制MPC要部署到云端或第三方高性能计算平台但状态量、控制量甚至模型参数都是敏感工业数据传统容错MPC在设计时又默认所有信息在控制器内部都是明文可见的。为了让连续搅拌式反应器CSTR这类典型连续生产装置在接入不可信计算平台之后既能保护隐私又能在传感器或执行器出故障时自动重构控制策略我在Matlab里试着把容错模型预测控制与同态加密从算法层面做了一次融合实现。这篇文章把整个思路、模型处理、代码结构和踩过的坑都写出来给正在做信息物理系统安全、工业MPC云化部署或者同态加密应用研究的同学一个可以直接参考的路线。先说明一点同态加密不是给你提供一个“更快的MPC”而是让MPC能在无法看到明文数据的机器上完成计算。这个前提决定了后面很多设计都要改变——包括控制器形式、求解器形式甚至故障检测模块放在哪一层。整条链路我在Matlab里全部跑通了一遍下面按研究动机、模型建立、加密求解改造、代码实现和性能权衡的顺序展开。1. 为什么非要把容错MPC和同态加密放在一起研究1.1 云化部署带来的真实隐私痛点模型预测控制在工业过程控制里的地位不用多说它能把约束处理、多变量耦合和未来趋势预测放在一个优化框架里解决。近年来不少工厂开始把MPC计算任务从本地DCS/PLC迁移到边缘服务器或云端原因很实际云端算力强可以跑更大规模、更短控制周期的优化多个产线可以共用一套控制器维护升级也集中。但问题随之而来。MPC每一拍都要读入当前状态$x_k$通常是温度、浓度、压力这类直接反映工艺运行情况的数据。这些数据在化工、制药、食品等行业直接关联配方和工艺Know-how很多企业不会允许它们离开生产网。更棘手的是MPC内部还有预测模型矩阵和约束矩阵这些参数一旦被第三方获取等于把装置的动态特性完整暴露了。传统做法是拉专线、签保密协议、做私有化部署但在多云协同、设备接入工业互联网平台的大趋势下总是有一些物理隔离做不到的场景。同态加密在这里提供了一个很有意思的思路如果计算平台永远接触不到明文状态只能看到加密后的密文同时还能在密文上执行加法、乘法运算让结果解密后等于明文计算结果那么MPC的核心优化计算理论上就可以外包给一个完全不可信的环境。这就是我在这项研究里选择同态加密作为隐私保护手段的根本原因——它不是防御式地“加密传输链路”而是让计算本身就在密文域完成。1.2 容错机制的前提竟然依赖明文状态容错MPC本身不是新概念。它做的事情是当传感器或执行器发生故障时整个闭环系统不能崩溃而是通过故障检测与诊断模块识别故障类型然后由控制器重构预测模型、约束条件或代价函数维持系统稳定并尽可能保住控制性能。典型主动容错MPC的闭环长这样标称MPC控制器给出控制量 → 执行器动作 → 传感器反馈量测 → 卡尔曼滤波/观测器生成残差 → 故障检测模块判断残差是否超限 → 如果确认故障重构模块修改MPC的某些参数。这里面有个默认前提卡尔曼滤波器和残差生成器需要明文状态MPC重构也需要明确的故障信息。换句话说传统容错MPC是建立在“控制器内部信息完全透明”基础上的。一旦引入同态加密尴尬就出现了。卡尔曼滤波里有矩阵乘法和求逆这些在密文域很难高效实现残差比较是典型的if-else逻辑同态加密下做比较非常昂贵MPC切换不同的故障模式本质上是一个决策树密文域很难支持这种分支计算。所以不能天真地把整个容错MPC原封不动搬进密文域必须重新划分安全边界哪些模块留在本地明文处理哪些模块可以加密外包这是整个研究里最核心的设计决策。1.3 融合设计的核心思想把计算逻辑“加密友好化”我在这项研究中最终采用了这样的分工本地保留状态估计、残差生成和故障决策因为这些模块强依赖比较、分支和解方程加密起来性价比极低云端或第三方平台只承接MPC优化求解的那部分计算且求解器必须改造成“不含比较、不含分支、只由加减乘除构成”的迭代形式。这样本地发出的不再是最原始的状态$x_k$而是$x_k$加密后的密文$\mathcal{E}(x_k)$云端回的也不直接是控制量$u_k$的明文而是若干轮迭代后的密文中间值本地解密后整理成最终控制量。这听起来好像只是把求解器换了个地方跑但实际工作量很大。因为商用MPC求解器大多基于内点法或有效集法每一步都要排序、比较、分支跳转同态加密下这些操作要么没法做要么代价高到不可接受。为了让MPC“适应”加密环境我不得不把优化问题改写成梯度下降迭代。后面第3章会详细讲这一步怎么改这里先记住结论同态加密环境下控制器必须是“迭代友好的”不是所有MPC形式都能直接被加密托管的。2. 连续搅拌式反应器的LTI模型与容错MPC落地表述2.1 在工作点附近把CSTR收紧成线性时不变系统连续搅拌式反应器本身是强非线性系统物料平衡和热量平衡里都带着Arrhenius指数项。但在工程上做MPC设计时通常不会直接在全局非线性模型上硬算而是在一个选定的稳态工作点附近做线性化得到一个线性时不变状态空间模型再围绕这个模型设计预测控制器。典型的CSTR两状态模型可以写成物料平衡$V\frac{dC_A}{dt} F(C_{Af}-C_A) - V k_0 e^{-E/(RT)} C_A$热量平衡$V\rho C_p\frac{dT}{dt} F\rho C_p(T_f-T) (-\Delta H)V k_0 e^{-E/(RT)} C_A - UA(T-T_c)$把状态选为$\Delta C_A$和$\Delta T$相对工作点的偏差输入选为冷却剂温度偏差$\Delta T_c$输出选为反应温度偏差$\Delta T$在工作点$(C_{As}, T_s, T_{cs})$附近做泰勒展开就能得到熟悉的LTI状态空间形式$$\dot{x} A x B u B_d d, \qquad y C x$$我在Matlab里用linmod或手工求Jacobian都做过实际效果差别不大。得到连续时间模型后再用c2d(sys, Ts)离散化Ts我取了0.1分钟6秒对应CSTR这种中等动态的化工对象是合理的。离散化之后的状态空间模型就是MPC的预测模型基础。实际做过CSTR控制的人都知道工作点选取直接影响线性化模型的精度。我建议不要选在开环不稳定点附近太近的位置否则后续加密域梯度下降迭代很容易因为模型本身条件数太大而发散。我的经验是把工作点选在反应转化率比较合适的稳态点并且在做完线性化后检查一下A矩阵特征值确保离散化后系统是可控可观的。2.2 故障场景定义与MPC重构策略容错MPC设计的第一步不是写代码而是明确要对付哪些故障。在这项研究中我重点考虑了三种工程上最常见的故障模式故障类型典型表现MPC重构策略执行器增益退化冷却阀实际动作幅度小于指令值修改输入矩阵B将B乘以估计的执行器效率因子执行器卡死冷却剂温度调节阀卡在某个开度约束中固定该输入剩余输入重新分配传感器偏置温度传感器读数偏离真实值在状态估计中补偿偏置或切换到冗余量测这三种故障都可以在LTI模型框架下表达。执行器增益退化最直接把原来的$B$矩阵替换为$\hat{B} B \cdot \alpha$$\alpha$是实时辨识出的效率因子执行器卡死则可以把对应输入通道从优化变量中剔除原先的输入约束变成一个固定值等式传感器偏置则在卡尔曼滤波的观测方程里加上一个可估计的偏置项。故障检测模块我用的是卡尔曼滤波器残差加滑动窗检验。具体做法用标称模型跑卡尔曼滤波得到新息序列残差$r_k y_k - C\hat{x}_{k|k-1}$正常情况下$r_k$是零均值白噪声故障发生时残差均值或方差会出现偏移。用一个长度为20的滑动窗计算残差均值再和阈值比较。这个阈值设多大会直接影响容错及时性和误报率我最后是通过蒙特卡洛仿真标定的后面避坑章节会说细节。MPC重构不需要像某些自适应控制那样在线辨识完整模型我的实现里是预设了几组故障模式对应的模型参数和约束集合故障检测模块触发后直接切换到相应配置。这样做的好处是故障切换是确定性的不会出现因为辨识不准导致控制器模型一直跳变的问题坏处是你得在离线阶段把主要故障场景想全这本身就是容错控制系统设计里绕不开的工程判断。2.3 MPC目标函数与约束的离散化表述在设计容错MPC时我用的目标函数是标准二次型$$J \sum_{j1}^{N_p} \hat{x}{kj|k}^T Q \hat{x}{kj|k} \sum_{j0}^{N_c-1} u_{kj}^T R u_{kj}$$预测时域$N_p$取20控制时域$N_c$取5。$Q$和$R$的选择比较常规$Q$对角元素根据温度和浓度的量级差做归一化$R$取一个相对较小的值让控制器不要太保守。约束方面考虑了执行机构物理限制冷却剂温度偏差$u \in [-2, 2]$控制增量$\Delta u \in [-0.5, 0.5]$。在每次采样周期MPC都需要解一个带约束的二次规划问题。这是整个融合研究里最关键也最麻烦的部分——因为带约束的QP求解器几乎都离不开比较和分支而这两个操作在密文域里恰恰是最贵的。为了让后面加密改造可行我没有直接调用quadprog而是先把约束做投影处理再用梯度投影法迭代求解。这样改造之后同一套MPC在明文环境跑和在加密环境跑求解器主体结构是基本一致的。3. 同态加密域中的MPC求解可行但不能照搬明文算法3.1 同态加密方案选型Paillier和CKKS到底能做什么同态加密并不是只有一种算法。做这个研究前我首先把方案选型想清楚了。目前最常接触的几类包括部分同态PHE如Paillier只支持加法同态ElGamal只支持乘法同态。近似同态SHE如BFV/BGV支持加法和乘法但电路深度有限。全同态FHE如CKKS支持浮点近似运算的加法和乘法可以做多层乘法电路。MPC优化求解需要什么运算只需要加法和乘法。理论上CKKS就够了。Paillier只能做加法那怎么处理乘法可以借助交互协议服务器把密文和一部分中间值发给客户端客户端解密后乘上另一个数再加密返回。这就是“外包计算交互”的思路。但交互次数多了通信开销就上去了。方案支持运算密文膨胀率计算速度适用场景Paillier加法同态约2048/明文位数较慢模幂运算需要交互式辅助乘法BFV/BGV整数加法和乘法约20-50倍中等整数编码场景CKKS浮点加法和乘法约10-20倍较快向量化浮点运算MPC很合适我在项目中先用Paillier把整个流程跑通因为它实现简单Matlab里手写也能跑适合理解整个交互链路后续为了提升计算效率又把核心求解部分切换到了CKKS思路不过Matlab生态对CKKS支持有限我是用C封装SEAL库再通过mex调用的。这里给个建议如果你只是验证算法可行性Paillier足够如果你想做性能评估建议直接上CKKS别在Paillier上浪费时间调优化。3.2 把带约束QP改造成“加密友好”的梯度投影迭代求解QP的办法有很多但能进密文域的很少。内点法要算牛顿方向、要做线搜索、要判断是否收敛这些操作在密文域里几乎没法高效实现。有效集法更依赖每一步的判断分支同样不适合。我在实现里采用的是梯度投影法核心迭代只有四步计算梯度$g 2\Phi^T Q (y_{ref} - \Phi u) 2 R u$这里的$\Phi$是预测模型构造出的分块矩阵。沿负梯度方向更新$u \leftarrow u - \eta g$。对$u$做投影把超出约束区间的分量拉回边界$u \min(\max(u, u_{lb}), u_{ub})$但密文域里min/max也不好做。这个问题的解法是把投影操作放回本地每个采样周期交互若干次由本地明文完成投影后重新加密上传。迭代若干次后停止取当前$u$作为控制量。第3步是重点。它意味着整个优化求解并不是“一锤子买卖”全部在云端完成而是在云端和本地之间做若干个回合的协作计算。这种模式在隐私计算里叫“server-aided computing”它不是纯非交互的同态计算安全性分析时要把交互信息一并考虑。好在这类交互暴露的只是投影后的中间值并且中间值每一轮都经过加密理论上不会泄露状态和约束敏感信息。收敛性问题也值得说。梯度投影法在普通约束QP上收敛速度中等但在加密域里不可能做自适应步长或者精确线搜索只能固定步长。步长设太大可能振荡设太小收敛很慢。我在实验中把步长$\eta$设为0.05迭代次数定为80次最终控制效果和quadprog相比性能损失在可接受范围内。3.3 浮点状态怎么变成整数并能完成运算同态加密不管是Paillier还是BFV在底层对浮点数的支持都很有限CKKS虽然支持浮点近似但众多意外噪声积累会导致精度下降。因此无论选哪个方案编码和量化都是绕不开的步骤。我的做法是这样的把状态变量先做归一化让取值范围落在$[-1,1]$区间。乘以一个缩放系数$s$并四舍五入取整。比如$s10^4$得到整数明文。对这个整数明文做同态加密。整个流程里缩放系数的选择直接影响精度和溢出风险。$s$太小量化误差大控制器输出噪声大$s$太大整数变大模运算溢出或密文规模膨胀的概率也增大。我建议对每个状态变量分别做归一化而不是统一缩放否则不同量纲温度偏差几度、浓度偏差零点几混在一起量化误差会互相放大。还有一个容易被忽略的问题是加减法后的缩放因子统一。状态乘一个系数、梯度公式里又除以一个系数如果缩放因子不一致明文编码下的数学关系和浮点运算就不是严格一致的了。我的做法是把整个MPC优化问题里的所有变量和系数都映射到同一个缩放空间里在明文调试阶段就固定好每个常数项的缩放倍数加密阶段不再改动。4. Matlab端到端实现方案与关键环节4.1 整体仿真框架与模块划分我在Matlab里的代码架构分成五个模块模型模块、故障注入模块、估计与检测模块、MPC求解模块、同态加密模块。整个仿真流程是模型模块定义CSTR的连续模型和离散模型用于模拟被控对象。故障注入模块在仿真过程中按预设时间点注入执行器增益退化、卡死或传感器偏置。估计与检测模块运行卡尔曼滤波器产生残差并做阈值判断。MPC求解模块接收故障检测结果选择对应的预测模型和约束运行梯度投影迭代。同态加密模块把MPC求解时需要的状态密文和中间结果做加密、解密和交互。这么做的好处是你可以在明文模式下先把容错MPC调好再切换到密文模式验证。如果一上来就加密出现问题根本分不清是控制算法的问题还是加密环节的问题。4.2 CSTR建模与离散化代码骨架CSTR模型的Matlab代码核心部分大概长这样% CSTR 线性化模型参数 A [ -1.341, -0.1689; -0.1789, -0.2087 ]; B [ 0; 0.3 ]; % 冷却剂温度输入 C [ 0, 1 ]; % 测量反应温度 D 0; Ts 0.1; % 采样周期 0.1 min sysC ss(A, B, C, D); sysD c2d(sysC, Ts, zoh); [Ad, Bd, Cd, Dd] ssdata(sysD);这里要说明一下实际A和B矩阵来自非线性CSTR方程在工作点的Jacobian我这里给的是简化参数只是为了说明代码结构。把注意力放在c2d的离散化方法上——我用了零阶保持器zoh因为是采样控制系统执行器在两个采样点之间保持恒定输出物理上更符合实际情况。MPC预测模型的构造用的是标准的延迟状态空间扩展。把状态增广成$[\Delta x; u_{k-1}]$就能把控制增量$\Delta u$当作新的优化变量。这样目标函数和约束都可以写成关于$\Delta u$的二次规划形式。这个扩展你可以在Matlab里用nargin之类的循环或直接分块矩阵生成本质就是一个分块Toeplitz矩阵的构造问题。4.3 卡尔曼滤波残差与故障重构的接入方式卡尔曼滤波器我用的是稳态增益版本省去在线递推Riccati方程的计算负担% 稳态卡尔曼增益 [~, L, ~] dare(Ad, Cd, Q_kf, R_kf); L L; % 这就是卡尔曼增益 K % 在每个采样周期 xp Ad * x_est Bd * u_prev; % 预测 res y_meas - Cd * xp; % 新息/残差 x_est xp L * res; % 更新残差$res$会进到一个滑动窗缓冲窗口长度取20个采样点。我对窗口内残差求均值计算一个简单统计量。超过阈值时就认为检测到了故障。阈值的选择我前面提到用蒙特卡洛标定在无故障情况下跑50次仿真统计残差均值的3倍标准差把阈值设为这个值。这样既保证误报率低又不至于让故障检测滞后太多。故障重构的接入点是MPC求解前的参数装配。我维护了一个结构体数组fault_config每一类故障对应一组B矩阵、约束上下限和代价权重。检测到故障后只需要把当前配置切换到对应条目下一拍MPC求解就自动用新模型。这个设计非常轻量也方便之后扩展新的故障类型。4.4 同态加密库在Matlab里怎么真正落地Matlab本身没有原生同态加密工具箱也没有官方支持的库。我试过两条路线第一条是手写Paillier加密。Paillier的核心就是模幂运算Matlab里用powermod能做但2048位模数下每个采样周期几十次模幂慢得无法接受。所以我只用它来验证正确性——确认加密、解密、加法的流程没有错确认量化编码的缩放因子一致。第二条是mex封装C库。我在Windows上编译了Microsoft SEAL封装了密钥生成、加密、乘法、加法、重线性化这几个函数然后通过mex编译成Matlab可调用的接口。这个路线更复杂但跑起来效率比Matlab原生大几十倍。如果你不想碰C也可以考虑用Java库比如Paillier库通过MATLAB的Java接口调用但Java和MATLAB之间的类型转换比较费事。无论哪条路线都要注意一件事不要把加密库的接口写进MPC求解器的核心循环。我的做法是定义了一个抽象接口HomomorphicEngine里面只有encrypt、decrypt、add、mul、rerand这几个方法。明文模式用一个假的PlainEngine实现直接返回原值密文模式才切换到真正的SEALEngine。这样调试控制算法时完全不碰加密需要验证加密正确性时再切换过去。5. 实验观察到的性能与安全权衡5.1 故障场景下的闭环响应我在三组场景下分别测试了纯明文容错MPC和密文容错MPC场景一正常运行无故障。场景二第50个采样时刻注入执行器增益退化增益从1退到0.6。场景三第80个采样时刻注入传感器偏置温度测量值突增0.3。场景二的结果最有代表性。明文MPC在故障检测模块确认故障后约3个采样周期内完成模型切换反应温度偏差最大超调约0.8密文MPC因为梯度迭代只跑到80步控制量更新本身就带一点收敛残差最大超调约1.1但仍在约束范围内闭环系统保持稳定。这说明加密改造没有破坏容错控制的基本目标主要代价在瞬态性能上。场景一无故障的对比更有意思密文MPC的稳态控制精度比明文MPC略低原因是量化误差和有限次梯度迭代造成了固定偏差。但不是完全失控稳态偏差大约在原值的2%以内。对于CSTR这种允许一定稳态偏差的化工对象这个精度是可以接受的。5.2 同态加密的计算开销到底有多大我分别测了四种配置下单步控制周期的耗时明文MPC、Paillier方案、CKKS方案SEAL实现、以及“明文检测密文求解”混合方案。采样周期为6秒结果如下配置单步控制耗时通信量/步说明明文MPC约0.02秒约几十字节本地方案性能基线明文检测Paillier密文求解约3.5秒约2 MB能跑但挤占采样周期明文检测CKKS密文求解约0.8秒约600 KB可实时但需要C封装全明文容错MPC约0.15秒—对比参照结论很清楚Paillier方案在实时性上很勉强CKKS则把单步耗时压缩到了1秒以内。如果你的控制周期只有0.1秒那这个方案还需要继续优化但如果控制周期是几秒甚至更长很多化工回路都是分钟级CKKS方案已经完全具备实时运行条件。我这里测的密钥规模是2048位、多项式阶数8192隐私强度足够工程使用。5.3 采样周期、迭代步数和密钥规模的联动选择这三个参数不是独立的。我把它们之间的关系整理成一个简单的工程判断链条采样周期越长留给同态运算的时间窗口越大密文规模可以选更大、迭代步数可以更多迭代步数越多控制精度越接近明文QL密钥规模越大单步耗时越长通信量也越大。所以实际选择时应该先定控制周期再反向倒推迭代步数和密钥规模。我自己的常用组合是控制周期6秒、CKKS多项式模数8192、迭代80步、步长0.05。如果控制周期压缩到3秒我会把迭代降到50步密钥规模不变控制精度损失大约5%左右。如果控制周期只有1秒我建议放弃纯加密方案改走“加密调度明文控制”的混合路边比如只在状态切换或批次切换时做加密验证。6. 实战中容易踩的坑与我的应对建议6.1 加密噪声导致的控制量抖动问题让我最头疼的问题在初次密文仿真时出现控制量每拍都在小幅波动看起来像高频抖动而不是平滑变化。排查后定位到两个来源一是量化误差在梯度计算里被放大尤其是$Q$矩阵权重大的通道量化噪声直接进入了梯度二是CKKS本身的近似噪声在多次乘法后累积导致同一输入在不同迭代轮次里得到略有不同的梯度值。解决办法是把目标函数里状态项做归一化时不要用太大缩放系数保持量化间隔适中更重要的是在明文编码阶段把权重矩阵$Q$和$R$也缩放到同一个整数空间避免加密域里出现小整数乘以大权重造成位数膨胀。实测下来把缩放系数从$10^4$降到$10^3$控制量抖动明显减少稳态精度损失不到1%。6.2 故障检测阈值和加密更新频率的配合问题还有一次系统误报让我反思了故障检测模块和加密更新频率的配合。因为云端MPC是每6秒解密一次状态、执行一轮优化本地卡尔曼滤波却是每0.5秒更新一次。状态估计频率比控制频率高得多导致同一控制周期内的残差不是平稳的滑动窗统计量偏大容易出现假报警。我的做法是把卡尔曼滤波的更新频率和控制器的执行频率强制对齐或者退一步在残差滑动窗里做时间加权让靠近当前时刻的残差权重更高。这样故障检测保持了时间上的敏锐度又不会因为高频估计噪声触发误报。6.3 从仿真走向工程之前必须想清楚的三件事这个研究做到后面我发现真正的瓶颈已经不在算法而在工程落地。第一密钥管理怎么做。同态加密不是对称加密随便发个密钥就行它涉及私钥、公钥、重线性化密钥、伽罗瓦密钥等多套密钥。私钥一旦泄露整个隐私保护形同虚设。工业现场没有专门的密钥管理服务这个部分要提早设计。第二云端计算结果的验证问题。既然云端不可信你凭什么相信它返回的密文控制量是正确计算出来的目前我依赖的是“解密后做边界合理性校验”但这并不能防恶意的错误计算。要真正解决这个问题需要再加零知识证明或可验证计算这部分我还在摸索。第三故障检测本身的数据隐私。我的实现里把残差生成留在了本地明文侧这是安全边界的合理选择。但如果你有多个厂区的数据都想上传做联合分析那本地残差统计量本身也是敏感信息需要单独设计隐私保护协议。最后一个个人体会收尾我最初做这个研究时也以为同态加密只是“接入一个库那么简单”实际走完一遍才发现真正的难点不是加密算法本身而是如何把一个成熟的、充满if-else和分支判断的控制算法改造成加密环境下能高效运行的形态。这个改造过程逼着我重新审视了MPC求解的每一步运算最终倒是把控制算法本身理解得更透彻了。如果你想复现这条路建议按“明文MPC → 明文迭代式MPC → 模拟密文接口 → 真加密实现”的顺序推进每一步都确认无误再进下一步。这样即便中途某个环节出了问题也能清楚定位在哪一层。