1. 线性MPC设计基础
任务目标
为一阶惯性系统设计线性MPC控制器，实现设定值跟踪。通过仿真对比不同预测时域 \(N_p\) 与控制时域 \(N_c\) 下的闭环响应，观察其对跟踪速度与超调量的影响。
系统模型与离散化
被控对象为一阶惯性系统（连续传递函数）：
\(G(s)=\dfrac{K}{Ts+1}\)，取 \(K=1,\ T=2\)
采样时间 \(T_s=0.5\text{s}\)，采用零阶保持器离散化，得到离散状态空间模型：
\(x(k+1)=a x(k)+b u(k)\)，\(y(k)=c x(k)\)
其中 \(a=e^{-T_s/T}=e^{-0.25}\approx0.7788\)，\(b=K(1-a)=0.2212\)，\(c=1\)。
MPC控制器设计
采用增量形式MPC，以控制增量 \(\Delta u\) 为决策变量，避免稳态误差。定义预测模型：
\(x(k+1)=a x(k)+b\bigl(u(k-1)+\Delta u(k)\bigr)\)
在预测时域 \(N_p\) 内迭代展开，构造预测输出向量 \(Y\) 与未来控制增量向量 \(\Delta U\) 的线性关系：
\(Y=F x(k)+G_u u(k-1)+G_{\Delta u}\Delta U\)
优化目标（设定值 \(r\) 跟踪，控制增量惩罚）：
\(\min_{\Delta U}\ (R-Y)^T Q (R-Y)+\Delta U^T R_w \Delta U\)
约束：控制增量幅值 \(|\Delta u|\le 0.5\)，控制量 \(u\in[-2,2]\)。该问题为二次规划（QP），调用 scipy.optimize.minimize 求解。
完整实现代码
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
# 系统参数
K, T, Ts = 1.0, 2.0, 0.5
a = np.exp(-Ts / T)
b = K * (1 - a)
c = 1.0
# MPC参数
Np = 10 # 预测时域
Nc = 3 # 控制时域
Q = 1.0 # 输出权重
Rw = 0.1 # 控制增量权重
du_max = 0.5 # 增量约束
u_max, u_min = 2.0, -2.0
# 构建预测矩阵
def build_matrices(Np, Nc, a, b, c):
# 状态预测: x(k+i) = a^i x(k) + sum_{j=0}^{i-1} a^{i-1-j} b u(k+j)
# 输出 y = c*x
# 构造 Y = F*x0 + G*U, 其中 U 为未来控制量序列
F = np.array([c * a**(i+1) for i in range(Np)])
G = np.zeros((Np, Np))
for i in range(Np):
for j in range(i+1):
G[i, j] = c * a**(i-j) * b
# 转换为增量形式: U = u(k-1)*ones + L*ΔU
L = np.tril(np.ones((Np, Np)))
G_du = G @ L[:, :Nc] # 只取前Nc个增量
G_u = G @ np.ones(Np) # u(k-1)的系数
return F, G_du, G_u
F, G_du, G_u = build_matrices(Np, Nc, a, b, c)
def mpc_controller(xk, u_prev, r, Np, Nc, Q, Rw, du_max, u_max, u_min):
# 预测输出: Y = F*xk + G_u*u_prev + G_du*ΔU
Y_free = F * xk + G_u * u_prev # 自由响应（ΔU=0）
def cost(dU):
Y = Y_free + G_du @ dU
e = r - Y
return e @ (Q * e) + Rw * dU @ dU
# 约束: |Δu|<=du_max, u_min<=u_prev+cumsum(ΔU)<=u_max
cons = []
# 增量约束
cons.append({'type': 'ineq', 'fun': lambda dU: du_max - np.abs(dU)})
# 控制量约束（累积）
def u_constraint(dU):
U_seq = u_prev + np.cumsum(dU)
return np.concatenate([U_seq - u_min, u_max - U_seq])
cons.append({'type': 'ineq', 'fun': u_constraint})
dU0 = np.zeros(Nc)
res = minimize(cost, dU0, method='SLSQP', constraints=cons,
options={'ftol': 1e-8, 'maxiter': 200})
return res.x[0] # 只取第一个增量
# 仿真
def simulate(Np, Nc, r=1.0, steps=60):
x = 0.0
u_prev = 0.0
y_hist, u_hist = [], []
for _ in range(steps):
du = mpc_controller(x, u_prev, r, Np, Nc, Q, Rw, du_max, u_max, u_min)
u = u_prev + du
u = np.clip(u, u_min, u_max)
x = a * x + b * u
y = c * x
y_hist.append(y)
u_hist.append(u)
u_prev = u
return y_hist, u_hist
# 对比不同时域
configs = [(5, 2), (10, 3), (20, 5)]
plt.figure(figsize=(10, 4))
for Np, Nc in configs:
F, G_du, G_u = build_matrices(Np, Nc, a, b, c)
y_hist, _ = simulate(Np, Nc)
plt.plot(y_hist, label=f'Np={Np}, Nc={Nc}')
plt.axhline(1.0, color='gray', ls='--', lw=0.8)
plt.xlabel('采样步数 k')
plt.ylabel('输出 y')
plt.legend()
plt.title('不同预测/控制时域下的设定值跟踪')
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('mpc_comparison.png', dpi=100)
plt.show()运行与验证
- 依赖安装：
pip install numpy scipy matplotlib - 运行：直接执行上述脚本，将生成
mpc_comparison.png并显示对比曲线。 - 通过条件：三条曲线均能从初始值0渐近收敛到设定值1.0；\(N_p=5,N_c=2\) 时响应最慢或出现小幅振荡，\(N_p=20,N_c=5\) 时跟踪最快但控制增量可能更频繁触及约束边界。若曲线发散或无法收敛，说明QP求解失败或矩阵构建有误。
时域影响分析
- 预测时域 \(N_p\)：增大 \(N_p\) 使控制器“看得更远”，稳定性增强，但过大会增加计算量并可能使响应过于保守。
- 控制时域 \(N_c\)：增大 \(N_c\) 提供更多自由度，跟踪更快，但过大会导致控制增量剧烈变化，系统易振荡。
实际应用中需根据工艺约束权衡二者：快速性要求高时适当增大 \(N_c\)，稳定性优先时增大 \(N_p\) 并限制 \(N_c\)。