4. MPC处理输出约束的软约束方法
任务目标
当MPC的输出约束过于严格或存在噪声时，硬约束可能导致优化问题无解。本案例采用软约束方法：允许输出在必要时轻微越界，但通过惩罚项抑制越界程度。设计一个单输入单输出（SISO）系统的MPC控制器，并分析约束惩罚权重 \( \lambda \) 对控制性能的影响。
系统模型与控制器设计
被控对象采用一阶惯性加纯延迟模型（FOPDT）：
\( G(s) = \frac{1}{10s+1} e^{-2s} \)
采样时间 \( T_s = 1\,\text{s} \)，离散化后得到状态空间模型（见代码中 sysd）。MPC采用增量型预测模型，预测时域 \( N_p = 10 \)，控制时域 \( N_c = 3 \)。输出约束为 \( y_{\min} = -0.5 \)，\( y_{\max} = 0.5 \)。软约束通过引入松弛变量 \( \epsilon \) 实现，优化目标为：
\[ \min_{\Delta u, \epsilon} \sum_{i=1}^{N_p} \| y_{k+i|k} - r \|^2_Q + \sum_{j=0}^{N_c-1} \| \Delta u_{k+j} \|^2_R + \lambda \epsilon^2 \]
满足：\( y_{\min} - \epsilon \le y_{k+i|k} \le y_{\max} + \epsilon \)，\( \epsilon \ge 0 \)。
关键完整代码
<code>import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from scipy.signal import lti, dlti, dstep
# 1. 系统模型（离散化）
Ts = 1.0
sys_c = lti([1], [10, 1]) # 一阶惯性
sys_d = dlti(sys_c, dt=Ts) # 零阶保持离散化
A_d, B_d, C_d, D_d = sys_d.to_ss().A[0,0], sys_d.to_ss().B[0,0], sys_d.to_ss().C[0,0], sys_d.to_ss().D[0,0]
# 2. MPC参数
Np, Nc = 10, 3
Q = 1.0
R = 0.1
lambda_pen = 50.0 # 软约束惩罚权重（可调）
y_min, y_max = -0.5, 0.5
u_min, u_max = -2.0, 2.0
du_max = 0.5
# 3. 预测矩阵构造（增量型）
# 状态: x = [y_k, u_{k-1}]^T, 输入增量 du
A = np.array([[A_d, B_d], [0, 1]])
B = np.array([[B_d], [1]])
C = np.array([1, 0])
# 构建预测矩阵 F, Phi
F = np.zeros((Np, 2))
Phi = np.zeros((Np, Nc))
for i in range(Np):
F[i] = C @ np.linalg.matrix_power(A, i+1)
for j in range(min(i+1, Nc)):
Phi[i, j] = C @ np.linalg.matrix_power(A, i-j) @ B
# 4. 软约束MPC优化函数
def mpc_soft(x0, r, u_prev):
# 决策变量: [du_0..du_{Nc-1}, epsilon]
n_du = Nc
def objective(z):
du = z[:n_du]
eps = z[n_du]
y_pred = F @ x0 + Phi @ du
cost = Q * np.sum((y_pred - r)**2) + R * np.sum(du**2) + lambda_pen * eps**2
return cost
# 约束: du限幅, epsilon>=0, 输出软约束
cons = []
# du限幅
cons.append({'type': 'ineq', 'fun': lambda z: du_max - np.abs(z[:n_du])})
# epsilon >= 0
cons.append({'type': 'ineq', 'fun': lambda z: z[n_du]})
# 输出软约束: y_min - eps <= y <= y_max + eps
def cons_out(z):
du = z[:n_du]
eps = z[n_du]
y_pred = F @ x0 + Phi @ du
return np.concatenate([y_pred - (y_min - eps), (y_max + eps) - y_pred])
cons.append({'type': 'ineq', 'fun': cons_out})
# 初始猜测
z0 = np.zeros(n_du + 1)
res = minimize(objective, z0, method='SLSQP', constraints=cons,
options={'ftol': 1e-8, 'maxiter': 200})
return res.x[:n_du], res.x[n_du]
# 5. 仿真
T_sim = 30
x = np.array([0.0, 0.0]) # [y, u_prev]
u_prev = 0.0
y_hist, u_hist, eps_hist = [], [], []
r_hist = np.ones(T_sim) * 0.8 # 设定值0.8，超出硬约束上限0.5
for k in range(T_sim):
r = r_hist[k]
du_opt, eps = mpc_soft(x, r, u_prev)
du = du_opt[0]
# 控制量限幅
u_new = np.clip(u_prev + du, u_min, u_max)
# 系统更新
y_new = A_d * x[0] + B_d * u_prev
x = np.array([y_new, u_new])
u_prev = u_new
y_hist.append(y_new)
u_hist.append(u_new)
eps_hist.append(eps)
# 6. 绘图
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(range(T_sim), y_hist, 'b-', label='y(k)')
plt.axhline(y=0.5, color='r', linestyle='--', label='y_max')
plt.axhline(y=-0.5, color='r', linestyle='--', label='y_min')
plt.plot(range(T_sim), r_hist, 'g--', label='r')
plt.ylabel('输出')
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(range(T_sim), u_hist, 'm-', label='u(k)')
plt.ylabel('控制量')
plt.xlabel('时间步 k')
plt.legend()
plt.tight_layout()
plt.savefig('soft_constraint_mpc.png', dpi=100)
print("仿真完成，结果图已保存为 soft_constraint_mpc.png")运行与验证
环境依赖：Python 3.8+，需安装 numpy、scipy、matplotlib。
运行方式：将上述代码保存为 soft_mpc.py，在终端执行 python soft_mpc.py。
验证方法：
- 观察输出曲线：设定值 \( r=0.8 \) 超出硬约束上限0.5，软约束允许输出短暂越界，但最终应稳定在0.5附近（因约束惩罚作用）。
- 检查
eps_hist：在初始阶段 \( \epsilon > 0 \)（允许越界），稳态后 \( \epsilon \to 0 \)（约束恢复）。 - 修改
lambda_pen值（如设为5和500），重新运行并对比：lambda_pen越大，输出越界幅度越小，但响应可能变慢；lambda_pen越小，越界越明显，但跟踪速度更快。
通过条件：程序无报错运行结束，生成 soft_constraint_mpc.png 图，且图中输出曲线在初始阶段出现轻微越界（超过0.5）后回落，控制量始终在 \([-2, 2]\) 内。
惩罚权重影响分析
软约束的核心是平衡约束满足与控制性能。\( \lambda \) 增大时，优化器更倾向于减小 \( \epsilon \)，即严格约束输出不越界，但可能牺牲跟踪速度或增加控制量波动；\( \lambda \) 减小时，系统允许较大越界以换取更快的设定值跟踪。实际应用中，\( \lambda \) 应根据工艺对约束违反的容忍度整定：安全关键约束取大值，非关键约束可适当减小。