💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。
⛳️座右铭：行百里者，半于九十。
📋📋📋本文目录如下：🎁🎁🎁
目录
💥1 概述
模型预测控制（MPC）的一个被广泛认可的缺点是，它通常只能在慢动态应用中使用，其中采样时间以秒或分钟为单位。实现快速MPC的一种众所周知的技术是离线计算整个控制律，在这种情况下，在线控制器可以作为查找表实现。此方法适用于状态和输入维度较小（例如，不超过五个）、约束较少且时间范围较短的系统。在本文中，我们描述了使用在线优化提高MPC速度的一系列方法。这些利用 MPC 问题的特定结构的自定义方法计算控制操作的速度比使用通用优化器的方法快 100 倍。例如，我们的方法在大约 12 毫秒内计算具有 3 个状态、30 个控制和 450 个时间步长（这需要求解具有 1284 个变量和 5 个约束的二次规划）的问题的控制动作，允许 MPC 在 200 Hz 下执行。
在经典模型预测控制（MPC）中，每个时间步的控制动作是通过求解在线优化问题获得的。对于线性模型、多面体约束和二次成本，生成的优化问题是二次规划 （QP）。使用通用方法求解QP可能很慢，这传统上将MPC限制在动态较慢的应用，采样时间以秒或分钟为单位。实现快速 MPC 的一种方法是将 QP 的解显式计算为初始状态的函数 [1]， [2];然后以查找表的形式在线实现控制操作。这里的主要缺点是表中的条目数量可以随着水平、状态和输入维度呈指数级增长，因此“显式 MPC”只能可靠地应用于小问题（其中状态维度不超过五个左右）。
在本文中，我们描述了一组方法，这些方法可用于使用在线优化大大加快MPC中控制动作的计算速度。一些想法已经在文献中被注意到，在这里我们将证明，当结合使用时，它们允许MPC的实现速度比通用优化器快几个数量级。
我们的主要策略是利用MPC [3]，[4]中出现的QP结构。已经注意到，通过适当的变量重新排序，可以通过求解线性方程组的块三对角线系统来找到每一步的内点搜索方向。利用这种特殊结构，状态维度的问题n、输入维度m和地平线T需要O(T(n+m)3)内点方法中每步骤的操作，而不是O(T3(n+m)3)如果不利用特殊结构。由于内点方法只需要恒定（且适度）的步骤数，因此MPC的复杂性在问题范围内是线性的，而不是立方的。
模型预测控制还有其他几个名称，例如滚动水平规划、后退水平控制、动态矩阵控制和动态线性规划。它已被广泛应用于化学过程和工业控制[5]，[9]-[12]，排队系统控制[13]，供应链管理[14]，[15]，经济和金融中的随机控制问题[16]，[17]，动态对冲[18]和收入管理[19]，[20]。
有关从最优控制角度讨论MPC，请参阅[21]–[23];有关稳定性和最优性结果的调查，请参见[24]。有关封闭式（脱机）方法，请参阅 [1]、[2] 和 [25]。有关使用优化方法的动态轨迹优化，请参见 [26]。
以前解决MPC中出现的QP有效解决方案的工作包括[27]-[29]和[63];有关化学过程控制中出现的大规模MPC问题的有效方法，请参阅[30]–[32]。[33]–[35]解决了稳健MPC的有效方法。使用牛顿方法应用于平滑非二次成本函数的想法可以在[36]，[37]中找到。
📚2 运行结果
部分代码：
%% Parameters
n = 8; % Dimension of state
m = 5; % Dimension of control
Q = eye(n); % State stage cost
R = eye(m); % Control stage cost
S = []; % State control coupled cost
Qf = 50*eye(n); % Terminal state cost
q = []; % Linear state cost
r = []; % Linear control cost
qf = []; % Terminal state cost
Xmax = 10; % State upper limit
Umax = 2; % Control upper limit
xmin = -Xmax*ones(n,1); % State lower bound
xmax = Xmax*ones(n,1); % State upper bound
umin = -Umax*ones(m,1); % Cotrol lower bound
umax = Umax*ones(m,1); % Control upper bound
high_limit = 1;
low_limit = 0;
A = (high_limit-low_limit).*rand(n,n) + ones(n,n)*low_limit; % Random A (State transition) matrix
B = (high_limit-low_limit).*rand(n,m) + ones(n,m)*low_limit; % Random B (Control matrix) matrix
A = A./(max(abs(eig(A)))); % Spectral radius of A within 1
high_limit_w = 1;
low_limit_w = 0;
w = (high_limit_w-low_limit_w).*rand(n,1) + ones(n,1)*low_limit_w; % Random noise vector
T = 10; % Horizon length
x0 = rand(n,1); % Initial state (random)
xf = 1*ones(n,1); % Terminal state
test = Fast_MPC2(Q,R,S,Qf,q,r,qf,xmin,xmax,umin,umax,T,x0,...
A,B,w,xf,[]); % Build class
%% Solve
% Native matlab solver
tic;
[x_opt_mat] = test.matlab_solve; % Solving using native matlab solver fmincon
t_mat = toc;
% Structured MPC full solve
tic;
[x_opt_full] = test.mpc_solve_full; % Solving structure problem as log barrier method with infeasible start newton
t_full = toc;
% Fixed log barrier method k=0.01
k_fix = 0.01;
tic;
[x_opt_log] = test.mpc_fixed_log(k_fix); % Fixed log(k) iteration method
figure(1);
stairs(x_mat);hold on;
stairs(x_full)
stairs(x_log);
stairs(x_nw);
stairs(x_lgnw);
legend('Matlab solver', 'In nw method', 'fixed log, fixed newton', 'fixed log + newton');
axis tight;
figure(2);
stairs(u_mat);hold on;
stairs(u_full)
stairs(u_log);
stairs(u_nw);
stairs(u_lgnw);
legend('Matlab solver', 'In nw method', 'fixed log, fixed newton', 'fixed log + newton');
axis tight;
%% Parameters
n = 8; % Dimension of state
m = 5; % Dimension of control
Q = eye(n); % State stage cost
R = eye(m); % Control stage cost
S = []; % State control coupled cost
Qf = 50*eye(n); % Terminal state cost
q = []; % Linear state cost
r = []; % Linear control cost
qf = []; % Terminal state cost
Xmax = 10; % State upper limit
Umax = 2; % Control upper limit
xmin = -Xmax*ones(n,1); % State lower bound
xmax = Xmax*ones(n,1); % State upper bound
umin = -Umax*ones(m,1); % Cotrol lower bound
umax = Umax*ones(m,1); % Control upper bound
high_limit = 1;
low_limit = 0;
A = (high_limit-low_limit).*rand(n,n) + ones(n,n)*low_limit; % Random A (State transition) matrix
B = (high_limit-low_limit).*rand(n,m) + ones(n,m)*low_limit; % Random B (Control matrix) matrix
A = A./(max(abs(eig(A)))); % Spectral radius of A within 1
high_limit_w = 1;
low_limit_w = 0;
w = (high_limit_w-low_limit_w).*rand(n,1) + ones(n,1)*low_limit_w; % Random noise vector
T = 10; % Horizon length
x0 = rand(n,1); % Initial state (random)
xf = 1*ones(n,1); % Terminal state
test = Fast_MPC2(Q,R,S,Qf,q,r,qf,xmin,xmax,umin,umax,T,x0,...
A,B,w,xf,[]); % Build class
%% Solve
% Native matlab solver
tic;
[x_opt_mat] = test.matlab_solve; % Solving using native matlab solver fmincon
t_mat = toc;
% Structured MPC full solve
tic;
[x_opt_full] = test.mpc_solve_full; % Solving structure problem as log barrier method with infeasible start newton
t_full = toc;
% Fixed log barrier method k=0.01
k_fix = 0.01;
tic;
[x_opt_log] = test.mpc_fixed_log(k_fix); % Fixed log(k) iteration method
figure(1);
stairs(x_mat);hold on;
stairs(x_full)
stairs(x_log);
stairs(x_nw);
stairs(x_lgnw);
legend('Matlab solver', 'In nw method', 'fixed log, fixed newton', 'fixed log + newton');
axis tight;
figure(2);
stairs(u_mat);hold on;
stairs(u_full)
stairs(u_log);
stairs(u_nw);
stairs(u_lgnw);
legend('Matlab solver', 'In nw method', 'fixed log, fixed newton', 'fixed log + newton');
axis tight;
🎉3 参考文献
文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。
566
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
