罗兰谈MATLAB的艺术

将想法转化为MATLAB

锥规划与最优离散动力学

发布的罗兰美国舒尔，

56次观看(过去30天)|喜欢|0评论
今天的客座博主是Alan Weiss，他为优化工具箱™和其他数学工具箱编写文档。

锥编程

嗨,各位。今天的主题是锥程序设计，以及锥程序设计在火箭优化控制中的应用。自R2020b以来 coneprog 求解器已用于求解锥规划问题。什么是锥程序设计?我认为它是二次规划的一种推广。所有二次规划问题都可以表示为锥规划问题。但是有些锥规划问题不能用二次规划来表示。
那么，什么是锥规划呢?这是一个具有线性目标函数和线性约束的问题，就像线性规划或二次规划。但是它也包含了锥约束。在三维[x, y, z]中，你可以将一个圆锥表示为，例如，一个圆在x-y方向上的半径小于或等于z。换句话说，圆锥约束就是不等式约束
$ x²+y²\le z²$ ，
或者同样的
$ \|[x,y]\|\le z $ 为非负 z ．
这是圆锥体边界的图片 $ \|[x,y]\|\le z $ 为非负 z ．
[X,Y] = meshgrid(-2:0.1:2);
Z =√(X。²+ y .²);
冲浪(X, Y, Z)
视图(8,2)
包含(“x”）
ylabel (“y”）
zlabel (“z”）
当然，您可以缩放、平移和旋转圆锥约束。一般圆锥约束的正式定义使用一个矩阵 Asc 向量, 二元同步通信 而且 d ，和标量 γ 约束条件为 x 表示为
norm(Asc*x - bsc) <= d'*x - gamma;
的 coneprog “优化工具箱”中的求解器要求您使用 secondordercone 函数来制定圆锥约束。例如,
Asc = diag([1,1/2,0]);
BSC = 0 (3,1);
D = [0;0;1];
Gamma = 0;
socConstraints = seconddordercone (Asc,bsc,d,gamma);
F = [- 1,2,0];
Aineq = [];
Bineq = [];
Aeq = [];
Beq = [];
lb = [-Inf，-Inf,0];
ub = [Inf,Inf,2];
[x,fval] = coneprog(f,socConstraints,Aineq,bineq,Aeq,beq,lb,ub)
找到最优解。
x = 3×1
0.4851 3.8806 2.0000
Fval = -8.2462
使用基于问题的方法来访问锥编程可能会更简单。该功能是在R2021a中添加的。对于前面使用基于问题的方法的例子:
X = optimvar(“x”，3,“下界”(负负0),“UpperBound”,正正2);
Asc = diag([1,1/2,0]);
问题=优化问题(“客观”x - x (1) 2 * (2));
概率。约束条件=范数(Asc*x) <= x(3);
[sol,fval] = solve(问题)
用coneprog解决问题。找到最优解。
索尔=带字段的结构:
X: [3×1 double]
Fval = -8.2462
请注意，与大多数非线性求解器不同，您不需要指定初始点 coneprog ．这在下面的例子中会派上用场。

圆锥约束的离散动力学

假设你想用最少的燃料控制火箭轻柔地降落在一个特定的位置。假设所使用的燃料与所施加的加速度乘以时间成正比。不要模拟燃烧燃料时火箭重量的变化;我们假设这个控制是相对较短的时间，权重没有明显的变化。在负z方向有重力加速度g = 9.81。火箭上也有线性阻力，在速度的负方向上，系数为1/10。这意味着时间过后 t ，在没有施加加速度或重力的情况下，速度从 v 来 $ v\exp(-t/10) $ ．
在连续时间中，位置的运动方程 p(t) $ 、速度 v(t) $ ，外加加速度 $ a(t) $ 是
$ \frac{dp}{dt} = v(t) $
$ \frac{dv}{dt} = -v(t)/10 + a(t) + g*[0,0，-1] $ ．
这里有一些近似的运动方程，使用离散时间 N 等距步长 $ t = t /N $ ：
$ p(i+i) = p(i) + t*(v(i) + v(i+1))/2 $ (梯形规则)
美元v (i + 1) =(我)* \ exp (- t / 10) + t *((我)+ g *[0, 0, 1])美元 (欧拉积分)。
因此,
$ p (i + 1) = p (i) + t * v (i) * (1 + \ exp (- t / 10)) / 2 + t ^ 2 *((我)+ g *[0, 0, 1]) / 2美元 ．
现在来看锥规划的部分。假设每一步所施加的加速度以常数为界 Amax ．这些约束条件是
$ \|a(i)\| \le {\rm Amax} $ 对所有 我 ．
最小化的代价应该是加速度的范数的和 t ．锥规划要求目标函数的优化参数为线性。您可以通过引入新的优化变量将此成本重新表述为线性的 (我) 它们受到一组新的锥约束:
$ {\rm cost} = \sum s(i)*t $
$ \|s(i)\| \le a(i) $ ．
假设火箭一开始以速度飞行 $ v0 = [100,50，-40] $ 在位置 $ p0 = [-1000，-800,1200] $ ．计算火箭到达预定位置所需的加速度 $ [0,0,0] $ 与速度 $ [0,0,0] $ 在时间 $ t = 40 $ ．将计算分为100个步骤( $ t=40/100 $ )．假设这是最大加速度 $ \rm{Amax} = 2g $ ．
的 makeprob 在 脚本结束 接受时间 T ，初始位置 p0 ，初速度 半 ，并返回一个描述离散动态和成本的问题。
P0 = [-1000，-800,1200];
V0 = [100,50，-40];
Prob = makprob (40,p0,v0)
概率=
带属性的优化问题:描述:" ObjectiveSense: '最小化'变量:[1×1 struct]包含4个优化变量目标:[1×1 OptimizationExpression]约束:[1×1 struct]包含4个优化约束参见问题公式与show。
设置选项以使用比默认值小100倍的最优容忍度来解决锥编程问题。使用 “舒尔” 线性求解器，可以更准确地解决这个问题。
Opts = optimoptions(“coneprog”，“OptimalityTolerance”1 e-8“LinearSolver”，“舒尔”）;
[sol,cost] = solve(问题，选项=opts)
用coneprog解决问题。找到最优解。
索尔=带字段的结构:
A: [99×3 double] p: [100×3 double] s: [99×1 double] v: [100×3 double]
成本= 312.7740
的 plottrajandaccel 函数 在这个脚本的末尾 画出轨道和加速度的范数作为时间步长的函数。
plottrajandaccel (sol)
最佳加速度几乎是“砰砰”。火箭在大约 $ 2g $ 开始时，加速度接近于零，直到接近终点。在接近终点时，火箭以最大速度加速以减缓下降速度并以零速度着陆。这个控制的总费用大约是313。

找到最佳时间

找出火箭着陆的最佳时间T，这意味着火箭使用尽可能少的燃料的时间。的 findT 函数 在这个脚本的末尾 调用 fminbnd 定位最小成本时间。我做了简单的实验，发现[20,60]是一个合理的时间范围 T 求最小值，我在 fminbnd 调用。如果你花的时间远远小于20，你会得到一个不可行的问题:
Badprob = makeprob(15,p0,v0);
badsol = solve(badprob,Options=opts)
用coneprog解决问题。问题是不可行的。
badsol =带字段的结构:
A: [] p: [] s: [] v: []
(题外话，如果你试着 T 一个优化变量，那么问题就不再是一个 coneprog 问题。相反，这是一个问题 fmincon ，在这种情况下需要更长的时间来解决，并且需要你提供一个初始点。)
Topt = findT(选项)
解决……完成
Topt = 22.3294
画出最佳轨迹和加速度。
probopt = makeprob(Topt,p0,v0);
[solopt,costopt] = solve(probopt,Options=opts)
用coneprog解决问题。找到最优解。
solopt =带字段的结构:
A: [99×3 double] p: [100×3 double] s: [99×1 double] v: [100×3 double]
成本值= 171.1601
plottrajandaccel (solopt)
最优成本约为171，大约是原始参数成本的一半。这一次，控制更接近bang-bang。火箭一开始加速到最大，然后停止加速一段时间。同样，在最后时刻，火箭以最大速度加速，以零速度着陆。

最终的想法

锥规划是解决许多凸优化问题的一个令人惊讶的通用框架。有关另一个重要示例，请参见 基于问题的圆锥规划最小化分段线性质量-弹簧系统能量 ．关于其他可以放在圆锥编程框架中的问题，请参阅Lobo, Miguel Sousa, Lieven Vandenberghe, Stephen Boyd和Hervé Lebret。二阶锥规划的应用 线性代数及其应用 284年,没有。1-3(1998年11月):193-228。 https://doi.org/10.1016/s0024 - 3795 (98) 10032 - 0
你觉得圆锥编程或离散动力学有用吗?你有什么自己的例子可以分享吗?让我们知道 在这里 ．

辅助函数

此代码创建 makeprob 函数。
函数轨迹问题= makeprob(T,p0,v0)
N = 100;
G = 9.81;
pF = [0 0 0];
Amax = 2*g;
P = optimvar(“p”N 3);
V = optimvar(“v”N 3);
A = optimvar(“一个”n - 1 3);
S = optimvar(“s”n - 1,“下界”0,“UpperBound”, Amax);
轨迹问题=优化问题;
t = t /N;
trajectoryproblem。目标= sum(s)*t;
scons = optimconstr(N-1);
为i = 1:(N-1)
Scons (i) = norm(a(i，:)) <= s(i);
结束
acons = optimconstr(N-1);
为i = 1:(N-1)
acons(i) = norm(a(i，:)) <= Amax;
结束
vcons = optimconstr(N+1,3);
Vcons (1，:) = v(1，:) == v0;
vcons (2: N:) = v (2: N:) = = v (1: (N - 1):) * exp (- t / 10) + t * (a + repmat ([0 0 - g), N - 1, - 1));
vcons(N+1，:) = v(N，:) == [0 0 0];
pcons = optimconstr(N+1,3);
Pcons (1，:) = p(1，:) == p0;
pcons (2: N:) = p (2: N:) = = p (1: (N - 1):) + (1 + exp (- t / 10)) / 2 * t * v (1: (N - 1):) + t ^ 2/2 * (a + repmat ([0 0 - g), N - 1, - 1));
pcons((N+1)，:) = p(N，:) == pF;
轨迹问题。约束。acons = acons;
trajectoryproblem.Constraints.scons = scons;
trajectoryproblem.Constraints.vcons = vcons;
轨迹问题。约束。pcons = pcons;
结束
此代码创建 plottrajandaccel 函数。
函数plottrajandaccel (sol)
数字
Psol = sol.p;
P0 = psol(1，:);
pF = psol(end，:);
plot3 (psol (: 1) psol (:, 2), psol (:, 3),“处方”）
持有在
plot3 (p0 (1), p0 (2), p0 (3),“ks”）
pF plot3 (pF (1), (2), pF (3),“波”）
持有从
视图([-10]18)
包含(“x”）
ylabel (“y”）
zlabel (“z”）
传奇(“步骤”，“起点”，“最后一点”）
数字
Asolm = sol.a;
Nasolm =根号(sum(asolm.^2,2));
情节(nasolm“处方”）
包含(“时间步”）
ylabel (“规范(加速度)”）
结束
此代码创建 fvalT 函数，由 findT ．
函数Fval = fvalT(T,opts)
P0 = [-1000，-800,1200];
V0 = [100,50，-40];
tprob = makeprob(T,p0,v0);
Opts = optimoptions(Opts，“显示”，“关闭”）;
[~，Fval] = solve(tprob,Options=opts);
结束
此代码创建 findT 函数。
函数Tmin = findT(opts)
disp (“解决……”）
Tmin = fminbnd(@(T)fvalT(T,opts)，20,60);
disp (“完成”）
结束
The MathWorks, Inc.版权所有
|
类别:
新特性,
优化
  • 打印
  • 发送电子邮件

另请参阅

评论

如欲留言，请点击在这里登录您的MathWorks帐户或创建一个新帐户。