【DR_CAN-MPC学习笔记】34.详细的MPC建模例子和matlab代码
上一篇博客:【DR_CAN-MPC学习笔记】2.最优化数学建模推导
参照二次规划一般形式,详细推导了MPC的数学模型,即最小化代价函数的表达式,最终推导结果为:
DR_CAN的视频:
【MPC模型预测控制器】3_一个详细的建模例子
【MPC模型预测控制器】3
【MPC模型预测控制器】4_完整案例讲解 - Octave代码
【MPC模型预测控制器】4
离散系统状态空间一般形式:
其中 为状态向量(n×1), 为输入向量(p×1), 为系统状态矩阵(n×n), 为系统输入矩阵(n×p)。
单输入二阶系统的例子:
上式中,,,n = 2,p = 1
系统状态向量和输入向量的关系:
- 表示在k时刻预测 k+1 时刻的系统状态。
- 由于 由 决定,因此不需要 ,所以 比 少一个维度。
- 因为初始值 和 均为 n×1 向量,因此 为 (N+1)n×1 向量。同理可推出 为 Np×1 向量。
- 为 (N+1)n×n 矩阵。
- 矩阵上面所有的 0 与初始状态 有关(n×1矩阵), 均为 n×p 矩阵,因此 为 (N+1)n×Np 矩阵。
具体可参考上一篇博客的推导(【DR_CAN-MPC学习笔记】2.最优化数学建模推导):
分析过程:
回到单输入二阶系统的例子,,,n = 2,p = 1,假设预测区间 N=3 ,
整理一下维度:
= | + | |||||
(N+1)n×1 | = | (N+1)n×n | n×1 | + | (N+1)n×Np | Np×1 |
8×1 | = | 8×2 | 2×1 | + | 8×3 | 3×1 |
对于系统输出方程: ,参考值 ,误差 ,代价函数为:
代价函数 = 误差加权和 + 输入加权和 + 终端误差,其中 和 为权重系数矩阵且均为对角矩阵。
经过简化后消去变量 (简化过程参考【DR_CAN-MPC学习笔记】2.最优化数学建模推导):
由上式可见, 只包含了初始状态项 和输入项 ,对 进行最优化可以得到输入项 。
矩阵 与 有关:
其中 和 是原来两个权重矩阵 和 的增广形式:
矩阵计算较为复杂,可用编程求解。
例子代码:
控制目标:设计合适的 使得 随着 的增加,趋近于0。引入误差 :
为目标值,控制目标为令误差接近0。
状态矩阵(n×k):
...
输入矩阵(p×k):
...
状态方程:
例如:k=1 时, 表示第2列
回到单输入二阶系统的例子:
为系统状态变量权重矩阵, 为系统输入变量权重矩阵, 为终端权重矩阵。
DR_CAN给出了Octave代码,在Matlab中也可以运行。下面的代码是我在此基础上修改后的单输入例子的代码,后面也有介绍如何修改为多输入系统。
传送门:(二输入系统)【MPC模型预测控制器】4_Octave代码
代码一共由三个部分组成,分别为主程序:MPC_Test.m,以及两个函数:MPC_Matrices.m和Prediction.m
MPC_Test.m
设置初始参数:
,...
%% 清屏
clear;
close all;
clc;%% 加载 optim package,若使用matlab,则注释掉此行
pkg load optim;%% 第一步,定义状态空间矩阵
%% 定义状态矩阵 A, n x n 矩阵
A = [1 0.1; 0 2];
n= size (A,1); % 计算矩阵第一个维度的长度%% 定义输入矩阵 B, n x p 矩阵
B = [0; 0.5];
p = size(B,2); % 计算矩阵第二个维度的长度%% 定义Q矩阵,n x n 矩阵
Q=[1 0; 0 1];%% 定义F矩阵,n x n 矩阵
F=[1 0; 0 1];%% 定义R矩阵,p x p 矩阵
R=[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]; % 初始状态不为0,控制目标为0%% 定义输入矩阵 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") %% 作者:DR_CAN https://www.bilibili.com/read/cv16891782 出处:bilibili
分析:
注释掉以下行,即输入设为0,可单独运行:
N=5;
[E,H]=MPC_Matrices(A,B,Q,R,F,N);
U_K(:,k) = Prediction(X_K(:,k),E,H,N,p);
得到:
由上图可见, 均趋于无穷,因为初始状态不为0且输入为0.
接下来设置合适的输入使得状态值趋于0.
代价函数:
与初始状态有关,不影响代价函数,因此控制目标为最小化 .
预测了N项,但结果只取第一项 .
矩阵 在之前已经分析过了:
矩阵的计算用 MPC_Matrices 函数解决,即代码:
[E,H]=MPC_Matrices(A,B,Q,R,F,N);
功能:输入矩阵 A,B,Q,R,F 和预测区间 N ,输出矩阵 E,H。具体过程参考下面的 MPC_Matrices.m 文件。接下来进行预测:
U_K(:,k) = Prediction(X_K(:,k),E,H,N,p);
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 矩阵, 得到 pM=[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 NP 个 0% 定义M 和 C
tmp=eye(n); %定义一个n x n 的 I 矩阵% 更新M和C
for i=1:N % 循环,i 从 1到 Nrows =i*n+(1:n); %定义当前行数,从i x n开始,共n行 C(rows,:)=[tmp*B,C(rows-n, 1:end-p)]; %将c矩阵填满tmp= A*tmp; %每一次将tmp左乘一次AM(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 %%作者:DR_CAN https://www.bilibili.com/read/cv16891782 出处:bilibili
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 = U_k(1:p,1); % 取第一个结果end %%作者:DR_CAN https://www.bilibili.com/read/cv16891782 出处:bilibili
分析:
quadprog:matlab自带的最优化函数
运行结果:
由上图所示,状态值 趋于0.
以上为单输入系统的例子。
二输入例子:
代码修改一下,也可以用于多输入,比如:
修改矩阵 A,B,R,Q,F,R ,使得 变化得更快(通过对矩阵Q的设置),且降低能耗减小 初始值(系统的输入一般是耗能的部分,通过对矩阵R的设置)
修改后的 MPC_Test.m:(其余不变)
%% 清屏
clear;
close all;
clc;%% 加载 optim package,若使用matlab,则注释掉此行
% pkg load optim;%% 第一步,定义状态空间矩阵
%% 定义状态矩阵 A, n x n 矩阵
% A = [1 0.1; 0 2];
A = [1 0.1; -1 2];
n= size (A,1); % 计算矩阵维度%% 定义输入矩阵 B, n x p 矩阵
% B = [0; 0.5];
B=[0.2 1; 0.5 2];
p = size(B,2);%% 定义Q矩阵,n x n 矩阵
% Q=[1 0; 0 1];
Q=[100 0; 0 1]; % 更加看重x_1的变化%% 定义F矩阵,n x n 矩阵
% F=[1 0; 0 1];
F=[100 0; 0 1];%% 定义R矩阵,p x p 矩阵
% R=[0.1];
R=[1 0; 0 0.1]; % 减小能耗,减小输入u_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") %% 作者:DR_CAN https://www.bilibili.com/read/cv16891782 出处:bilibili
由上图所示, 迅速趋近于0. 相比于 R=[0.1 0; 0 0.1] 时,R=[1 0; 0 0.1] 时的输入 减小,但 趋近于0的速度变缓。由此也可以看出MPC的最优化的决策结果不是绝对的。
PS:
腾讯云开学季优惠:
【腾讯云】多款云产品1折起,买云服务器送免费机器,最长免费续3个月
【腾讯云】爆款云服务器限时体验20元起,更多上云必备产品低至1元
【DR_CAN-MPC学习笔记】34.详细的MPC建模例子和matlab代码相关推荐
- JUC.Condition学习笔记[附详细源码解析]
JUC.Condition学习笔记[附详细源码解析] 目录 Condition的概念 大体实现流程 I.初始化状态 II.await()操作 III.signal()操作 3个主要方法 Conditi ...
- 观察者模式学习笔记(详细)
观察者模式学习笔记(详细) 一.什么是观察者模式 观察者模式,是定义对象之间的一对多的关系,主要作用是减少对象之间的耦合度,分为两个角色 被观察者:其实就是发布者,发布消息通知所有的观察者 观察者:接 ...
- 吴恩达《机器学习》学习笔记四——单变量线性回归(梯度下降法)代码
吴恩达<机器学习>学习笔记四--单变量线性回归(梯度下降法)代码 一.问题介绍 二.解决过程及代码讲解 三.函数解释 1. pandas.read_csv()函数 2. DataFrame ...
- 【Unity学习笔记】[Unity中文课堂教程] C#中级编程代码
[Unity学习笔记][Unity中文课堂教程] C#中级编程代码 最近想补一补C#基础,Unity官方的C#中级编程教程质量很高,于是开个帖子把跟着敲+记录了部分价讲解和我自己的理解的代码存在这 原 ...
- MPC学习笔记1:基于状态空间模型的预测控制(1)
MPC调节器 1.给定一个由状态空间法描述的离散系统: MPC控制器与其他线性二次调节器(LQR)的区别就在于其可以很好的将系统动态约束纳入考虑. 采样周期Ts控制了算法的效率,太大会错过很多系统运行 ...
- MPC学习笔记(1)——原理
最近在学习M. W. Mehrez的MPC时发现了很多不了解的细节,分享一下对该算法的梳理与理解. 在自动驾驶或机器人领域中,模型预测控制(Model Predictive Control, MPC) ...
- SpringBoot学习笔记整理详细
写在前面:欢迎来到「发奋的小张」的博客.我是小张,一名普通的在校大学生.在学习之余,用博客来记录我学习过程中的点点滴滴,也希望我的博客能够更给同样热爱学习热爱技术的你们带来收获!希望大家多多关照,我们 ...
- 影像组学视频学习笔记(34)-使用3D Slicer软件提取影像组学特征、Li‘s have a solution and plan.
作者:北欧森林 链接:https://www.jianshu.com/p/afcd06221ea4 来源:简书,已获授权转载 RadiomicsWorld.com "影像组学世界" ...
- DR_CAN的学习笔记--1现代控制理论
本文是在学习B站up主DR_CAN博士之后的学习笔记以及一些个人的感悟.B站链接为https://space.bilibili.com/230105574/,感兴趣的可以去B站看视频. 本文针对的读者 ...
最新文章
- windows多线程同步--临界区
- golang append时slice len 和 cap
- 正则表达式验证代码(字母、数字、Email、网址、电话号码、汉字、身份证号码)
- Python实训day14am【Python网络爬虫综合大作业-答辩】
- 实验-网页动画(js版)
- PWN-PRACTICE-BUUCTF-17
- 智能机器人语音识别技术详细解析
- SVN下载及语言包安装
- matlab对矩阵谱分解
- maya加载不了arnold的mtoa可能是这个低级错误!
- 量化金融分析AQF(1):股票概述
- Olivetti PR2/PR2E 打印机故障分析与排除
- CPU时钟周期和时钟频率
- 用计算机画小鸡,水墨电脑画--丝瓜小鸡图
- java常用二进制数据转换工具
- html与css第三天
- 关于编程C++——如何写程序
- java二进制补码_java基础 二进制补码
- 用云服务器windows环境下来搭建一个Minecraft服务器教程以及客户端使用教程java版
- Win+ubuntu单硬盘双系统安装(雷神新911机型)