基于高斯(Gaussian)过程的数据驱动模型预测控制算法GP-MPC 包括：高斯(Gaussian)过程GP代码；多状态高斯(Gaussian)过程GP模型预测控制MPC代码；单状态高斯(Ga

985计算机硕士

于 2025-05-04 09:45:00 发布

阅读量853

点赞数 28

文章标签： python django pygame virtualenv flask tornado scikit-learn

本文链接：https://blog.csdn.net/2501_90824910/article/details/147620059

版权

基于高斯(Gaussian)过程的数据驱动模型预测控制算法GP-MPC

包括：高斯(Gaussian)过程GP代码；
多状态高斯(Gaussian)过程GP模型预测控制MPC代码；
单状态高斯(Gaussian)过程GP模型预测控制MPC代码；
实例1线性系统

文章目录

高斯过程（Gaussian Process, GP）是一种非参数化的贝叶斯方法，广泛应用于回归、分类以及作为数据驱动的模型预测控制（MPC）中的动态模型。GP 能够通过训练数据学习复杂的非线性系统动力学，并提供具有不确定性的预测，这对于 MPC 的鲁棒性和安全性非常有用。

下面我将分别展示：

高斯过程（GP）的基本代码实现（使用 scikit-learn）
单状态和多状态系统的高斯过程模型预测控制（GP-MPC）
实例 1：线性系统 - 直流电机模型
实例 2：非线性系统 - Van der Pol 振荡器

一、高斯过程（GP）回归基础代码

我们首先演示一个基本的 GP 回归示例，使用 scikit-learn 库。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C

# 生成样本数据
def f(x):
    return np.sin(2 * np.pi * x)

X_train = np.random.rand(20, 1) * 2 - 1  # 训练输入 [-1, 1]
y_train = f(X_train).ravel() + 0.1 * np.random.randn(20)  # 带噪声输出

# 定义核函数
kernel = C(1.0, (1e-3, 1e3)) * RBF(1.0, (1e-2, 1e2))

# 创建 GP 回归模型
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10)
gp.fit(X_train, y_train)

# 预测
X_test = np.linspace(-1, 1, 100).reshape(-1, 1)
y_pred, sigma = gp.predict(X_test, return_std=True)

# 绘图
plt.figure()
plt.plot(X_test, y_pred, 'b-', label="Prediction")
plt.fill_between(X_test.ravel(), y_pred - 1.96 * sigma, y_pred + 1.96 * sigma,
                 color='gray', alpha=0.3, label="95% Confidence Interval")
plt.scatter(X_train, y_train, c='r', label="Training Data")
plt.legend()
plt.title("Gaussian Process Regression")
plt.show()

二、单状态系统的 GP-MPC 控制器

假设系统为一个单变量动态系统：

$x_{k+1} = f(x_k, u_k)$

我们使用 GP 学习该映射 $f$ ，然后将其嵌入到 MPC 中进行优化。

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor

class GPMPC_SingleState:
    def __init__(self, gp_model, N=10, umax=1.0, umin=-1.0):
        self.N = N  # 预测步数
        self.gp = gp_model
        self.umax = umax
        self.umin = umin

    def predict(self, x, u_seq):
        """给定初始状态和控制序列，预测未来的状态序列"""
        X = []
        x_pred = x.copy()
        for u in u_seq:
            x_pred = self.gp.predict(np.array([x_pred, u]).reshape(1, -1))[0]
            X.append(x_pred)
        return np.array(X)

    def cost(self, u_seq, x0, target):
        """MPC 成本函数"""
        X = self.predict(x0, u_seq)
        cost = np.sum((X - target)**2)  # 状态误差平方
        cost += 0.1 * np.sum((u_seq)**2)  # 控制量惩罚
        return cost

    def optimize(self, x0, target):
        """优化控制序列"""
        from scipy.optimize import minimize
        u_guess = np.zeros(self.N)
        res = minimize(lambda u: self.cost(u, x0, target), u_guess,
                       bounds=[(self.umin, self.umax)]*self.N)
        return res.x[0]  # 返回第一个最优控制输入

# 示例用法：
# 先训练 GP 模型（略），然后初始化 GPMPC_SingleState

三、多状态系统的 GP-MPC 控制器

对于多状态系统：

$\mathbf{x}_{k+1} = \mathbf{f}(\mathbf{x}_k, \mathbf{u}_k)$

我们使用多个 GP 分别建模每个状态维度，或使用多输出 GP。

from sklearn.multioutput import MultiOutputRegressor

# 示例：多状态 GP 模型
class MultiStateGP:
    def __init__(self):
        self.model = MultiOutputRegressor(GaussianProcessRegressor())

    def fit(self, X, Y):
        self.model.fit(X, Y)

    def predict(self, X_in):
        return self.model.predict(X_in)

# 多状态 GP-MPC
class GPMPC_MultiState:
    def __init__(self, gp_model, N=10, umax=1.0, umin=-1.0):
        self.N = N
        self.gp = gp_model
        self.umax = umax
        self.umin = umin

    def simulate(self, x0, u_seq):
        X = [x0.copy()]
        for u in u_seq:
            xu = np.hstack((X[-1], u))
            next_state = self.gp.predict(xu.reshape(1, -1))[0]
            X.append(next_state)
        return np.array(X)

    def cost_func(self, u_seq, x0, target):
        X = self.simulate(x0, u_seq)
        state_cost = np.sum((X[:-1] - target)**2, axis=1).sum()
        control_cost = 0.1 * np.sum(u_seq**2)
        return state_cost + control_cost

    def optimize_control(self, x0, target):
        from scipy.optimize import minimize
        u0 = np.zeros(self.N)
        res = minimize(lambda u: self.cost_func(u, x0, target), u0,
                       bounds=[(self.umin, self.umax)]*self.N)
        return res.x[0]

四、实例 1：直流电机 GP-MPC

4.1 系统模型

直流电机的动态方程如下：

$\dot{\omega} = \frac{K_t}{J} i - \frac{b}{J} \omega \\ \dot{i} = \frac{1}{L}(v - R i - K_e \omega)$

其中：

$\omega$ : 角速度
$i$ : 电流
$v$ : 输入电压
参数设为默认值

4.2 使用 GP-MPC 进行控制

# 生成训练数据
def dc_motor_dynamics(t, state, v, params):
    omega, i = state
    Kt, Ke, J, b, L, R = params
    domega = (Kt / J) * i - (b / J) * omega
    di = (1 / L) * (v - R*i - Ke*omega)
    return [domega, di]

params = [0.01, 0.01, 0.1, 0.1, 0.5, 1.0]  # Kt, Ke, J, b, L, R
from scipy.integrate import solve_ivp

# 生成训练数据
num_samples = 100
X_data = []
Y_data = []

for _ in range(num_samples):
    state = np.random.uniform([-10, -1], [10, 1])  # omega, current
    v = np.random.uniform(-12, 12)  # voltage input
    sol = solve_ivp(lambda t, s: dc_motor_dynamics(t, s, v, params), [0, 0.1], state, dense_output=True)
    next_state = sol.sol(0.1)
    X_data.append(np.hstack((state, v)))
    Y_data.append(next_state)

X_train = np.array(X_data)
Y_train = np.array(Y_data)

# 训练 GP 模型
multi_gp = MultiStateGP()
multi_gp.fit(X_train, Y_train)

# 初始化 MPC
mpc = GPMPC_MultiState(multi_gp, N=10, umax=12, umin=-12)

五、实例 2：Van der Pol 振荡器 GP-MPC

5.1 Van der Pol 动态模型

$\ddot{x} - \mu(1 - x^2)\dot{x} + x = u$

5.2 GP-MPC 实现

def van_der_pol(t, state, u, mu=1.0):
    x, dx = state
    ddx = mu*(1 - x**2)*dx - x + u
    return [dx, ddx]

# 生成训练数据
X_vdp = []
Y_vdp = []
for _ in range(100):
    state = np.random.uniform([-2, -2], [2, 2])
    u = np.random.uniform(-1, 1)
    sol = solve_ivp(lambda t, s: van_der_pol(t, s, u), [0, 0.1], state, dense_output=True)
    next_state = sol.sol(0.1)
    X_vdp.append(np.hstack((state, u)))
    Y_vdp.append(next_state)

# 训练 GP
multi_gp_vdp = MultiStateGP()
multi_gp_vdp.fit(X_vdp, Y_vdp)

# 初始化 MPC
mpc_vdp = GPMPC_MultiState(multi_gp_vdp, N=10, umax=1.0, umin=-1.0)

总结

✅ 提供了 单状态和多状态 GP 模型及其在 MPC 中的应用
✅ 包括两个系统示例：直流电机（线性） 和 Van der Pol 振荡器（非线性）
✅ 所有代码均可扩展为带噪声的环境，只需在训练数据中加入噪声即可

在这里插入图片描述
高斯过程模型预测控制（GP-MPC），并生成与您提供的图像相似的结果，我们需要完成以下几个步骤：

定义直流电机的动力学模型：用于生成训练数据和仿真。
使用高斯过程回归拟合动力学模型。
设计并实现基于GP的MPC控制器。
仿真并绘制结果。

以下是完整的Python代码示例，包括上述所有步骤：

import numpy as np
from scipy.integrate import solve_ivp
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C
import matplotlib.pyplot as plt

# 定义直流电机的动力学模型
def dc_motor_dynamics(t, state, u, params):
    omega, theta = state
    Kt, Ke, J, b, L, R = params
    domega = (Kt / J) * u - (b / J) * omega
    dtheta = omega
    return [domega, dtheta]

# 生成训练数据
def generate_training_data(params, num_samples=100):
    X_train = []
    Y_train = []
    for _ in range(num_samples):
        state = np.random.uniform([-10, -np.pi], [10, np.pi])
        u = np.random.uniform(-12, 12)
        sol = solve_ivp(lambda t, s: dc_motor_dynamics(t, s, u, params), [0, 0.1], state, dense_output=True)
        next_state = sol.sol(0.1)
        X_train.append(np.hstack((state, u)))
        Y_train.append(next_state - state)
    return np.array(X_train), np.array(Y_train)

# 训练高斯过程模型
def train_gp_model(X_train, Y_train):
    kernel = C(1.0, (1e-3, 1e3)) * RBF(1.0, (1e-2, 1e2))
    gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10)
    gp.fit(X_train, Y_train)
    return gp

# GP-MPC 控制器
class GPMPC:
    def __init__(self, gp_model, N=10, umax=12, umin=-12):
        self.N = N
        self.gp = gp_model
        self.umax = umax
        self.umin = umin

    def predict(self, x, u_seq):
        X = [x.copy()]
        for u in u_seq:
            delta_x = self.gp.predict(np.array([X[-1][0], X[-1][1], u]).reshape(1, -1))[0]
            next_state = X[-1] + delta_x
            X.append(next_state)
        return np.array(X)

    def cost(self, u_seq, x0, target):
        X = self.predict(x0, u_seq)
        state_cost = np.sum((X[:, 0] - target)**2)  # 只考虑角速度误差
        control_cost = 0.1 * np.sum(u_seq**2)
        return state_cost + control_cost

    def optimize(self, x0, target):
        from scipy.optimize import minimize
        u_guess = np.zeros(self.N)
        res = minimize(lambda u: self.cost(u, x0, target), u_guess,
                       bounds=[(self.umin, self.umax)]*self.N)
        return res.x[0]

# 参数设置
params = [0.01, 0.01, 0.1, 0.1, 0.5, 1.0]  # Kt, Ke, J, b, L, R

# 生成训练数据并训练GP模型
X_train, Y_train = generate_training_data(params)
gp_model = train_gp_model(X_train, Y_train)

# 初始化GP-MPC控制器
mpc = GPMPC(gp_model, N=10, umax=12, umin=-12)

# 仿真
target_omega = 1.0
num_steps = 35
x = np.array([0.0, 0.0])  # 初始状态：角速度和角度
u_gp_mpc = []
x_gp_mpc = [x]
for k in range(num_steps):
    u = mpc.optimize(x, target_omega)
    u_gp_mpc.append(u)
    sol = solve_ivp(lambda t, s: dc_motor_dynamics(t, s, u, params), [0, 0.1], x, dense_output=True)
    x = sol.sol(0.1)
    x_gp_mpc.append(x)

x_gp_mpc = np.array(x_gp_mpc)
u_gp_mpc = np.array(u_gp_mpc)

# 绘图
fig, axs = plt.subplots(3, 1, figsize=(10, 8))

axs[0].plot(range(num_steps+1), x_gp_mpc[:, 0], 'g-', label='GP-MPC')
axs[0].set_title('Angular Velocity')
axs[0].legend()

axs[1].plot(range(num_steps+1), x_gp_mpc[:, 1], 'g-', label='GP-MPC')
axs[1].set_title('Angle')
axs[1].legend()

axs[2].plot(range(num_steps), u_gp_mpc, 'b-', label='GP-MPC')
axs[2].set_title('Input')
axs[2].legend()

plt.tight_layout()
plt.show()