MPC模型预测控制原理和Matlab以及Python代码实现

MPC模型预测控制原理和代码

一. 介绍模型预测控制(MPC)原理

在这里插入图片描述

简要解释一下最优控制最优控制的目标是在一定的约束条件下达到最优的系统表现,那么要让系统达到最优表现,一般是通过定义损失函数J,通过最小化损失函数J来达到最优控制,对于单入单出(SISO)系统来说,损失函数J上面已经定义了,多入多出(MIMO)系统的损失函数和SISO系统的区别就是单入单出系统的损失函数里面的q和r是实数,MIMO系统的损失函数J里面的Q和R是矩阵,它们用来控制误差u和输入u的权重,e和u分别代表误差和输入;

在这里插入图片描述在这里插入图片描述在这里插入图片描述

第二和第三部分是推导部分,自己推一推吧,没法解释很清楚;

在这里插入图片描述

第四部分是状态空间方程的矩阵的维度的计算;

上面就是MPC模型预测控制的原理,主要还是要自己推导一下,下面给出模型预测控制的Python代码以及matlab代码;

二. Python代码

import numpy as np
import scipy.linalg
from cvxopt import solvers, matrix
import matplotlib.pyplot as plt

A = np.array([[1, 1], [-1, 2]])
n = A.shape[0]

B = np.array([[1, 1],[1, 2]])
p = B.shape[1]

Q = np.array([[1, 0],[0, 1]])
F = np.array([[1, 0],[0, 1]])
R = np.array([[1, 0],[0, 0.1]])

k_steps = 100

X_k = np.zeros((n, k_steps))

X_k[:,0] = [10, -10]

U_k = np.zeros((p, k_steps))

N = 5

def cal_matrices(A,B,Q,R,F,N):

    n = A.shape[0]
    p = B.shape[1]

    M = np.vstack((np.eye((n)), np.zeros((N*n,n))))
    C = np.zeros(((N+1)*n,N*p))
    tmp = np.eye(n)

    for i in range(N):
        rows = i * n + n
        C[rows:rows+n,:] = np.hstack((np.dot(tmp, B), C[rows-n:rows, 0:(N-1)*p]))
        tmp = np.dot(A, tmp)
        M[rows:rows+n,:] = tmp

    Q_bar_be = np.kron(np.eye(N), Q)
    Q_bar = scipy.linalg.block_diag(Q_bar_be, F)
    R_bar = np.kron(np.eye(N), R)

    G = np.matmul(np.matmul(M.transpose(),Q_bar),M)
    E = np.matmul(np.matmul(C.transpose(),Q_bar),M)
    H = np.matmul(np.matmul(C.transpose(),Q_bar),C) + R_bar

    return H, E

def Prediction(M,T):

    sol = solvers.qp(M,T)
    U_thk = np.array(sol["x"])
    u_k = U_thk[0:2, :]

    return u_k

M, C = cal_matrices(A,B,Q,R,F,N)
M = matrix(M)

for k in range(1,k_steps):
    x_kshort = X_k[:, k - 1].reshape(2, 1)
    u_kshort = U_k[:, k - 1].reshape(2, 1)
    T = np.dot(C,x_kshort)
    T = matrix(T)
    for i in range(2):
        U_k[i:,k-1] = Prediction(M,T)[i,0]

    X_knew = np.dot(A,x_kshort) + np.dot(B,u_kshort)

    for j in range(2):
        X_k[j:,k] = X_knew[j,0]

print(X_k)

三 .Matlab代码

%% 加载 optim package,若使用matlab,则注释掉此行

pkg load optim;

%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 第一步,定义状态空间矩阵

%% 定义状态矩阵 A, n x n 矩阵

A = [1 0.1; -1 2];

n= size (A,1);

%% 定义输入矩阵 B, n x p 矩阵

B = [ 0.2 1; 0.5 2];

p = size(B,2);

%% 定义Q矩阵,n x n 矩阵

Q=[100 0;0 1];

%% 定义F矩阵,n x n 矩阵

F=[100 0;0 1];

%% 定义R矩阵,p x p 矩阵

R=[1 0 ;0 .1];

%% 定义step数量k

k_steps=100; 

%% 定义矩阵 X_K, n x k 矩 阵

X_K = zeros(n,k_steps);

%% 初始状态变量值, n x 1 向量

X_K(:,1) =[20;-20];

%% 定义输入矩阵 U_K, p x k 矩阵

U_K=zeros(p,k_steps);


%% 定义预测区间K

N=5;

%% Call MPC_Matrices 函数 求得 E,H矩阵 

[E,H]=MPC_Matrices(A,B,Q,R,F,N);


%% 计算每一步的状态变量的值

for k = 1 : k_steps 

%% 求得U_K(:,k)

U_K(:,k) = Prediction(X_K(:,k),E,H,N,p);

%% 计算第k+1步时状态变量的值

X_K(:,k+1)=(A*X_K(:,k)+B*U_K(:,k));

end


%% 绘制状态变量和输入的变化

subplot  (2, 1, 1);

hold;

for i =1 :size (X_K,1)

plot (X_K(i,:));

end

legend("x1","x2")

hold off;

subplot (2, 1, 2);

hold;

for i =1 : size (U_K,1)

plot (U_K(i,:));

end

legend("u1","u2")

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

~~~~~~~~~~~~~~~~~~~~~~~MPC_Matrices.m~~~~~~~~~~~~~~~~~~~~

function  [E , H]=MPC_Matrices(A,B,Q,R,F,N)


n=size(A,1);   % A 是 n x n 矩阵, 得到 n

p=size(B,2);   % B 是 n x p 矩阵, 得到 p

%%%%%%%%%%%%

M=[eye(n);zeros(N*n,n)]; % 初始化 M 矩阵. M 矩阵是 (N+1)n x n的, 

                         % 它上面是 n x n 个 "I", 这一步先把下半部

                         % 分写成 0 

C=zeros((N+1)*n,N*p); % 初始化 C 矩阵, 这一步令它有 (N+1)n x NP0

% 定义MC 

tmp=eye(n);  %定义一个n x n 的 I 矩阵

% 更新M和C

for i=1:N % 循环,i 从 1N

    rows =i*n+(1:n); %定义当前行数,从i x n开始,共n行 

    C(rows,:)=[tmp*B,C(rows-n, 1:end-p)]; %将c矩阵填满

    tmp= A*tmp; %每一次将tmp左乘一次A

    M(rows,:)=tmp; %M矩阵写满

end 


% 定义Q_bar和R_bar

Q_bar = kron(eye(N),Q);

Q_bar = blkdiag(Q_bar,F);

R_bar = kron(eye(N),R); 


% 计算G, E, H

G=M'*Q_bar*M; % G: n x n

E=C'*Q_bar*M; % E: NP x n

H=C'*Q_bar*C+R_bar; % NP x NP 


end

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

~~~~~~~~~~~~~~~~~~~~~~~Prediction.m~~~~~~~~~~~~~~~~~~~~~

function u_k= Prediction(x_k,E,H,N,p)

U_k = zeros(N*p,1); % NP x 1

U_k = quadprog(H,E*x_k);

u_k = U_k(1:p,1); % 取第一个结果

end

这就是模型预测控制MPC的原理以及代码实现,欢迎大家观看学习。

  • 22
    点赞
  • 109
    收藏
    觉得还不错? 一键收藏
  • 8
    评论
【资源说明】 1、该资源包括项目的全部源码,下载可以直接使用! 2、本项目适合作为计算机、数学、电子信息等专业的课程设计、期末大作业和毕设项目,作为参考资料学习借鉴。 3、本资源作为“参考资料”如果需要实现其他功能,需要能看懂代码,并且热爱钻研,自行调试。 基于Matlab实现模型预测控制(MPC)源码+项目说明.zip <center>MPC模型预测控制</center> = [toc] # [资料1. DR_CAN:MPC模型预测控制器](https://www.bilibili.com/read/cv16891782) MPC_Test.m,Copy_of_MPC_Matrics.m,Copy_of_Prediction.m是DR_CAN提供的示例代码,例子中的参考信号R=0,且输出方程y=x,即矩阵C为单位阵。([代码地址](https://www.bilibili.com/read/cv16891782)) MPC_demo.mlx 以一个二阶系统为例演示MPC,对DR_CAN的代码进行了拓展,参考信号可设,输出方程中的c可设。([笔记推导](./MPC_notes.pdf)) 在控制的教材中,常常考虑参考信号为0的简化情况。在参考信号不为0的情况下,可以通过引入误差$e=z-z_{d}$,将误差作为新的状态量,可以将问题重新转换为参考信号=0的情况(参考误差信号为0),这在 资料3.无人驾驶车辆模型预测控制 的推导笔记中有所涉及,引入了误差,同时通过误差实现了非线性系统的线性化。 # [资料2. MATLAB中国【Model Predictive Control】](https://space.bilibili.com/1768836923/search/video?keyword=mpc) ## 1. 特点 [参考视频1](https://www.bilibili.com/video/BV16U4y1c7EG?spm_id_from=333.999.0.0&vd_source=be5bd51fafff7d21180e251563899e5e) [参考视频2](https://www.bilibili.com/video/BV1Qu411Z7DQ/?spm_id_from=333.788.recommend_more_video.-1&vd_source=be5bd51fafff7d21180e251563899e5e) ### 1.1. 优点 1. 可以处理MIMO,而PID只能处理SISO,虽然可以使用多个PID控制多个变量,但当变量之间存在耦合时,PID参数的调节会很困难; 2. 可以处理约束条件,由于模型预测控制是通过构建优化问题来求解控制器的动作的,所以可以非常自然的将这些约束建立在优化问题中以此来保证这些约束的满足。; 3. 使用了未来的预测信息。 ### 1.2. 缺点 要求强大的计算力,因为在每一个时间步都需要求解优化问题。 ## 2. 参数设置 采样时间:设置为开环系统响应上升时间Tr的1/20~1/10 $$\frac{Tr}{20} \leqslant Ts \leqslant \frac{Tr}{10}$$ 预测区间:20~30个时间步 控制区间:预测区间的10%~20%,并且至少有2~3个时间步 约束条件:约束分为硬约束和软约束,硬约束不可违背,软约束可以违背。不建议对输入和输出都进行硬约束,因为两者可能冲突以致无法求解优化问题。建议将输出设为软约束,并避免对输入和输入变化率都有硬约束 权重:取决于实际情况 [参考视频](https://www.bilibili.com/video/BV1b44y1v7Xt/?spm_id_from=333.788.recommend_more_video.-1&vd_source=be5bd51fafff7d21180e251563899e5e) ## 3. 自适应MPC,增益调度MPC,非线性MPC 适用于处理非线性系统,其中自适应MPC和增益调度MPC的本质是将系统线性化 [参考视频](https://www.bilibili.com/video/BV1ZL411g7Ya/?spm_id_from=333.788.recommend_more_video.-1&vd_source=be5bd51fafff7d21180e251563899e5e) ### 3.1. 自适应MPC(Adaptive MPC) 处理非线性系统时,在每个工作点附近对系统作线性化,得到一个新的线性模型,使用的前提是优化问题的结构在每个工作点不变,即在约束范围内,状态数量和约束数量不变。 ### 3.2. 增益调度MPC(Gain-

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 8
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值