matlab中的mpc模型代码

下述在matlab中编写的mpc模型预测控制代码
图注标出画线部分不明白代码含义

img

b站 DR_CAN 编写的代码
具体mpc文字笔记详见链接https://www.robotsfan.com/posts/fe8d7b17.html
请解答，谢谢

% 清屏
clear;
close all;
clc;
% 第一步，定义状态空间矩阵
% 定义状态矩阵 A, n*n 矩阵
A = [1 0.1; -1 2];
n = size(A,1);
% 定义输入矩阵 B, n*p 矩阵
B = [0.2 1;0.5 2];
p = size(B,2);
% 定义Q矩阵，n*n 矩阵
Q = [100 0;0 1];
% 定义F矩阵，n*n 矩阵
F = [100 0;0 1];
% 定义R矩阵，p*p 矩阵
R = [1 0;0 0.1];
% 定义step数量k
k_steps = 100;
% 定义矩阵 X_K， n*k 矩 阵
X_K = zeros(n,k_steps);
% 初始状态变量值， n*1 向量
X_K(:,1) = [20;-20];
% 定义输入矩阵 U_K， p*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

function [E,H] = MPC_Matrices(A,B,Q,R,F,N)
    n=size(A,1); % A是n*n矩阵,得到n
    p=size(B,2); % B是n*p矩阵,得到p
    M=[eye(n);zeros(N*n,n)]; % 初始化M矩阵,M矩阵是(N+1)n*n的,
                             % 它上面是n*n个"I",这一步先把下半部分写成0
    C=zeros((N+1)*n,N*p); % 初始化C矩阵,这一步令它有(N+1)n*NP个0
    % 定义M和C
    tmp=eye(n); % 定义一个n*n 的 I 矩阵
    % 更新Ｍ和C
    for i=1:N % 循环,i从1到N
        rows =i*n+(1:n); %定义当前行数,从i*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*n
    E=C'*Q_bar*M; % E: NP*n
    H=C'*Q_bar*C+R_bar; % NP*NP




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

该回答引用ChatGPT

for i=1:N % 循环,i从1到N
        rows =i*n+(1:n); %定义当前行数,从i*n开始，共n行
        C(rows,:)=[tmp*B,C(rows-n, 1:end-p)]; %将c矩阵填满
        tmp= A*tmp; %每一次将tmp左乘一次A
        M(rows,:)=tmp; %将M矩阵写满

这一段代码是 MPC_Matrices 函数中的 for 循环，用于更新控制矩阵 C 和状态矩阵 M。具体来说，循环中的代码含义如下：

1、对于第 i 个循环，计算出 C 矩阵中的第 rows 行到第 rows+n-1 行，即 rows = i*n+(1:n)。其中，i 表示当前的循环次数，n 是状态矩阵 A 的行数，即状态变量的数量。这一行代码的作用是将控制矩阵 C 中第 i 个时刻到第 i+N 个时刻的输入矩阵和之前的状态矩阵组合在一起，得到一个 Np 行的控制矩阵，其中 p 是输入矩阵 B 的列数，即控制输入的数量。

2、接着，将当前的 tmp 矩阵（初值为单位矩阵）左乘状态矩阵 A，得到新的 tmp 矩阵。这一行代码的作用是在状态矩阵 A 的作用下，更新当前时刻的状态变量对未来时刻的影响。

3、最后，将更新后的 tmp 矩阵写入 M 矩阵中的第 rows 行到第 rows+n-1 行，即 M(rows,:) = tmp。这一行代码的作用是将更新后的状态矩阵写入 M 矩阵中，得到一个 N+1 行的状态矩阵，其中第 1 行为当前时刻的状态矩阵，第 2 到第 N+1 行为未来时刻的状态矩阵。

该回答引用GPTᴼᴾᴱᴺᴬᴵ
这是一个基于状态空间模型的MPC控制器的MATLAB代码实现。其中，通过定义状态空间矩阵A、输入矩阵B、权重矩阵Q、目标矩阵F和约束矩阵R，通过MPC_Matrices函数计算出E、H矩阵。然后通过Prediction函数进行状态预测，并将预测值作为优化问题的输入，计算出控制器输出值。最后通过循环迭代，计算出每个时间步的状态变量和输入的变化，并将结果绘制成图形展示。

该MPC控制器可以应用于多种系统控制，通过调整输入矩阵B、权重矩阵Q和约束矩阵R等参数，可以实现不同系统的控制需求。

它C通过连接两个其他矩阵tmp*B和来填充一个矩阵C(rows-n, 1:end-p)。

该代码使用循环迭代i从 1 到 的索引N，其中N是某个正整数。对于每个i，它定义了一个行索引向量rows，该向量从 开始in+1并包括n之后的连续整数（即，in+1到i*n+n）。

该行用两个矩阵的串联填充矩阵中由索引C(rows,:)=[tmpB,C(rows-n, 1:end-p)];的行。第一个矩阵是和的乘积，表示为。第二个矩阵是比当前行 ( ) 更早开始行的子矩阵。该语法意味着包括该子矩阵的所有列，但最后一列除外。这两个矩阵一起水平连接并分配给由in索引的行。rowsCtmpBtmpBCnC(rows-n, 1:end-p)1:end-pprowsC

参考GPT和自己的思路，这部分代码是MPC_Matrices函数中的一部分，其目的是构建控制问题的约束矩阵C。下面对这两句代码进行解释：

rows = in+(1:n); % 定义当前行数，从in开始，共n行
C(rows,:) = [tmp*B,C(rows-n,1:end-p)]; % 将C矩阵填满

首先，由于预测区间为N，因此需要将状态和控制输入向量的维度进行扩展。C矩阵的每一行对应了一个状态向量和其对应的控制输入向量。因此，C矩阵的维度应该为(N+1)n x Np。

对于每一个时间步i，我们需要将对应的状态向量和控制输入向量添加到C矩阵中。在这个过程中，我们使用了一个临时变量tmp来存储状态矩阵A的幂次。具体而言，我们在第i次迭代中将tmp乘以A，并将结果作为下一次迭代的tmp。同时，我们将tmp乘以B并将结果添加到C矩阵中，同时将上一次迭代的C矩阵中的对应行添加到当前行中。这个过程可以使用下面的公式表示：

C(rows,:)=[tmp*B,C(rows-n, 1:end-p)];

这里，rows变量表示当前行的索引，由于每个状态向量和控制输入向量都包含n和p个元素，因此它们的总维度为n+p。因此，我们使用C(rows,:)来选择当前行，然后将两个向量水平连接起来。tmpB表示将tmp与B矩阵相乘得到的一个np的矩阵。而C(rows-n,1:end-p)表示C矩阵上一次迭代中对应的行，从第一列到第(end-p)列。这个过程可以保证C矩阵中的每一行都包含了当前时间步之前的所有状态和控制输入。