模型预测控制：线性MPC

模型预测控制（Model Predictive Control, MPC）是一种广泛应用于工业过程控制和自动驾驶等领域的先进控制技术。MPC通过在线解决优化问题来计算控制输入，从而实现系统的最优控制。本文将介绍线性MPC的系统模型、优化问题、LMPC算法，以及单输入系统和多输入系统的具体应用。

系统模型（离散时间线性时不变 (LTI) 系统）

在线性MPC中，系统模型通常假设为离散时间线性时不变（LTI）系统。LTI系统的状态空间模型可以表示为：

$x_{k+1} = Ax_k + Bu_k$

其中：

$x_k$ 是时间步 $k$ 的系统状态向量
$u_k$ 是时间步 $k$ 的控制输入向量
$A$ 是状态转移矩阵
$B$ 是输入矩阵

输出方程通常为：

$y_k = Cx_k + Du_k$

其中：

$y_k$ 是时间步 $k$ 的输出向量
$C$ 是输出矩阵
$D$ 是直接传输矩阵

优化问题

在MPC中，控制输入是通过解决一个在线优化问题来确定的。优化问题的目标是最小化一个代价函数（通常是一个二次型函数），并满足系统的约束条件。典型的代价函数包括以下两部分：

状态跟踪误差：衡量实际状态与目标状态之间的差距。
控制输入：控制输入的大小，避免过大的控制动作。

代价函数通常可以表示为：

$\sum_{i=0}^{N-1} \left( x_{k+i|k}^T Q x_{k+i|k} + u_{k+i|k}^T R u_{k+i|k} \right) + x_{k+N|k}^T P x_{k+N|k}$

其中：

$N$ 是预测时域长度
$Q$ 是状态误差权重矩阵
$R$ 是控制输入权重矩阵
$P$ 是终端状态权重矩阵

优化问题需要满足系统的动态方程和控制输入的约束：

$\begin{align*} x_{k+i+1|k} &= A x_{k+i|k} + B u_{k+i|k} \quad \forall i = 0, \ldots, N-1 \\ x_{min} &\leq x_{k+i|k} \leq x_{max} \quad \forall i = 1, \ldots, N \\ u_{min} &\leq u_{k+i|k} \leq u_{max} \quad \forall i = 0, \ldots, N-1 \end{align*}$

LMPC算法

线性模型预测控制（LMPC）算法的步骤如下：

状态估计：获取当前状态 $x_k$ 。
优化求解：解决当前时刻的优化问题，得到预测时域内的最优控制输入序列 $\{u_{k|k}, u_{k+1|k}, \ldots, u_{k+N-1|k}\}$ 。
应用控制：应用当前时刻的控制输入 $u_k = u_{k|k}$ 。
滚动优化：将预测时域向前推进一个时间步长，重复以上步骤。

单输入系统

对于单输入单输出（SISO）系统，模型和控制输入的表示会更加简单。假设系统状态向量 $x_k$ 为一维向量，控制输入 $u_k$ 也是一维标量。系统模型可以表示为：

$x_{k+1} = ax_k + bu_k$

输出方程为：

$y_k = cx_k + du_k$

优化问题与多输入系统相同，只是矩阵 $A$ 、 $B$ 、 $C$ 、 $D$ 变成了标量。

多输入系统

对于多输入多输出（MIMO）系统，状态向量 $x_k$ 、控制输入向量 $u_k$ 和输出向量 $y_k$ 均为多维向量。系统模型可以表示为：

$x_{k+1} = A x_k + B u_k$

输出方程为：

$y_k = C x_k + D u_k$

MIMO系统的优化问题和单输入系统类似，只是需要处理更高维度的矩阵。

代码示例

以下是一个简单的Python代码示例，展示如何实现单输入单输出系统的线性MPC：

import numpy as np
from scipy.linalg import solve_discrete_are
import matplotlib.pyplot as plt

# 定义系统模型
A = np.array([[1.0]])  # 状态转移矩阵
B = np.array([[1.0]])  # 控制输入矩阵
C = np.array([[1.0]])  # 输出矩阵
D = np.array([[0.0]])  # 直通项矩阵

# 定义MPC参数
Q = np.array([[1.0]])  # 状态权重矩阵
R = np.array([[1.0]])  # 控制权重矩阵
N = 10  # 预测时域长度

# 离散时间代数Riccati方程求解P矩阵
P = solve_discrete_are(A, B, Q, R)

# 计算K矩阵
K = np.linalg.inv(B.T @ P @ B + R) @ (B.T @ P @ A)

# 目标设置
x_target = np.array([[5]])  # 目标状态
u_target = 0  # 目标控制输入

# 初始化状态
x = np.array([[0.0]])  # 初始状态

# 仿真时间
T = 50
x_trajectory = []
u_trajectory = []

for t in range(T):
    # 计算控制输入
    u = -K @ (x - x_target) + u_target
    
    # 应用控制输入并更新状态
    x = A @ x + B @ u
    
    # 记录状态和控制输入
    x_trajectory.append(x.item())
    u_trajectory.append(u.item())

# 绘制结果
plt.figure(figsize=(12, 6))

plt.subplot(2, 1, 1)
plt.plot(x_trajectory, label='State x')
plt.axhline(y=x_target.item(), color='r', linestyle='--', label='Target x')
plt.title('State Trajectory')
plt.xlabel('Time step')
plt.ylabel('State')
plt.legend()

plt.subplot(2, 1, 2)
plt.plot(u_trajectory, label='Control Input u')
plt.axhline(y=u_target, color='r', linestyle='--', label='Target u')
plt.title('Control Input Trajectory')
plt.xlabel('Time step')
plt.ylabel('Control Input')
plt.legend()

plt.tight_layout()
plt.show()

在这里插入图片描述

代码说明：

系统模型定义：
- A、B、C、D 分别为状态转移矩阵、控制输入矩阵、输出矩阵和直通项矩阵。
- 对于SISO系统，所有矩阵都是1x1的。
MPC参数：
- Q 和 R 是权重矩阵，用于调整状态和控制输入在优化中的重要性。
- N 是预测时域的长度，决定了MPC预测未来多少步。
离散时间代数Riccati方程求解：
- 使用scipy.linalg.solve_discrete_are函数求解P矩阵。
- 计算反馈增益矩阵K。
目标设置：
- x_target 是目标状态。
- u_target 是目标控制输入。
初始化状态：
- x 是初始状态。
仿真循环：
- 在每个时间步t，计算控制输入u，更新状态x，并记录状态和控制输入的轨迹。
结果绘制：
- 使用matplotlib库绘制状态和控制输入的轨迹。