罗兰的MATLAB艺术

将想法转化为MATLAB

锥形编程和最佳离散动态

作者罗兰美国舒尔，

56 次查看（过去 30 天） | 个赞 | 0个评论
今天的客座博主是Alan Weiss，他为Optimization Toolbox™和其他数学工具箱编写文档。

锥形编程

嗨伙计。今天的主题是锥形编程，以及锥形编程的应用，以最佳地控制火箭。自R2020B以来 coneprog 求解器已经可以解决锥形编程问题。什么是锥形编程？我认为它是二次编程的概括。所有二次编程问题都可以表示为锥形编程问题。但是有锥形编程问题不能代表为二次程序。
什么是圆锥规划?它是一个具有线性目标函数和线性约束的问题，就像一个线性规划或二次规划。但它也包含锥约束。在三维[x, y, z]中，你可以将一个圆锥表示为，例如，一个圆在x-y方向上的半径小于或等于z。换句话说，圆锥约束就是不等式约束
$ x^2+y^2\le z^2 $ ，
或者同样的
$ \ | [x，y] \ | \ le z $ 为非负 z ．
这是一张圆锥边界的图片 $ \ | [x，y] \ | \ le z $ 为非负 z ．
(X, Y) = meshgrid (2:0.1:2);
Z =√X。^ 2 + y ^ 2);
冲浪(X, Y, Z)
查看（8,2）
Xlabel（“x”）
ylabel（“y”）
Zlabel（“z”）
当然，您可以扩展，翻译和旋转锥限制。一般锥形约束的正式定义使用矩阵 Asc ，vectors. 二元同步通信 和 d 和标量 γ 在约束条件下 x 表示为
norm(Asc*x - bsc) <= d'*x - gamma;
的 coneprog “优化工具箱”中的求解器要求您使用 secondordercone 功能制定锥限制。例如，
Asc =诊断接头([1,1/2,0]);
二元同步通信= 0 (3,1);
d = [0, 0, 1];
γ= 0;
socConstraints = secondordercone (Asc,二元同步通信,d,γ);
f = [-1，-2,0];
Aineq = [];
bineq = [];
Aeq = [];
说真的= [];
磅=(负负0);
乌兰巴托=(正,正无穷,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”，[inf，Inf，2]）;
Asc =诊断接头([1,1/2,0]);
prob = OptimProblem（“客观的”x - x (1) 2 * (2));
概率。约束条件=范数(Asc*x) <= x(3);
[溶胶，fval] =解决（prob）
使用coneprog解决问题。找到最优解。
sol =结构体字段:
x(3×1双):
fval = -8.2462
请注意，与大多数非线性求解器不同，您不需要指定初始点 coneprog ．这在下面的示例中很有用。

具有锥限制的离散动态

假设您希望使用最小燃料控制火箭在特定位置轻轻地降落。假设所使用的燃料与应用的加速时间成比例。在燃烧燃料时，不要模拟火箭的变化重量;我们认为，这种控制是相对较短的时间，其中重量不会明显变得明显。在负Z方向上有重力加速度G = 9.81。在火箭上也存在线性拖动，其在具有系数1/10的负速的负方向上起作用。这意味着经过时间 t ，在没有施加任何加速度或重力的情况下，速度从 v 来 v \ exp (- t / 10)美元 ．
在连续时间中，运动方程的位置 $ p（t）$ 、速度 v (t)美元 ，并施加加速度 美元(t)美元 是
$ \ frac {dp} {dt} = v（t）$
$ 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, -40);
概率= makeprob (p0, v0)
概率=
描述:" ObjectiveSense: '最小化'变量:[1×1 struct]包含4个OptimizationVariables目标:[1×1 OptimizationExpression]约束:[1×1 struct]包含4个OptimizationConstraints参见show问题公式。
设置选项以使用比默认值小的100倍来解决锥形编程问题。使用 “舒尔” 线性求解器，可以更准确的解决这个问题。
选择= optimoptions (“coneprog”，“OptimalityTolerance”1 e-8“LinearSolver”，“舒尔”）;
(溶胶、成本)=解决(概率,选择=选择)
使用coneprog解决问题。找到最优解。
sol =结构体字段:
答：[99×3双] P：[100×3双] S：[99×1双] V：[100×3双]
成本= 312.7740
的 plottrajandaccel 函数 在这个脚本的末尾 画出加速度的轨迹和范数作为时间步长的函数。
plottrajandaccel (sol)
最佳加速度几乎是“砰砰”。火箭加速到大约 2 g美元 起初，然后在近端直到零加速度接近。靠近末端，火箭最大加速，以减慢零速度的下降和落地。该控件的总成本约为313。

找到最佳的时间

找出火箭着陆的最佳时间T，即使火箭使用尽可能少燃料的时间。的 findT 函数 在这个脚本的末尾 调用 fminbnd 找到最小的成本时间。我简单地尝试发现[20,60]是次数的合理范围 T 最低，我使用了这些界限 fminbnd 称呼。如果您花了一段时间不到20，则获得一个不可行的问题：
badprob = makeprob (15 p0 v0);
badsol =解决(badprob选项=选择)
使用coneprog解决问题。问题是不可行的。
badsol =结构体字段:
A: [] p: [] s: [] v: []
(顺便说一句，如果你想 T 那么问题不再是一个优化变量 coneprog 问题。相反，它是一个问题 fmincon ，在本例中需要更长的时间来解决，并且需要您提供一个初始点。)
topt = findt（选择）
解决……完成
Topt = 22.3294
绘制最佳的轨迹和加速度。
probopt = makeprob (p0, Topt v0);
[solopt, costopt] =解决(probopt选项=选择)
使用coneprog解决问题。找到最优解。
solopt =结构体字段:
答：[99×3双] P：[100×3双] S：[99×1双] V：[100×3双]
costopt = 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 函数。
函数traometoryproblus = 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);
trajectoryproblem = optimproblem;
t = t / N;
trajectoryproblem。目标= (s)和* t;
您= 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;
trajectoryproblem.Constraints.acons = acons;
trajectoryproblem.Constraints.scons =您;
trajectoryproblem.constraints.vcons = Vcons;
trajectoryproblex.constraints.pcons = pcons;
结束
此代码创建 plottrajandaccel 函数。
函数plottrajandaccel (sol)
数字
psol = sol.p;
: p0 = psol (1);
pF = psol(结束:);
Plot3（PSOL（：，1），PSOL（:,2），PSOL（:,3），“处方”）
持有在
plot3 (p0 (1), p0 (2), p0 (3),“ks”）
Plot3（PF（1），PF（2），PF（3），'博'）
持有从
视图([-10]18)
Xlabel（“x”）
ylabel（“y”）
Zlabel（“z”）
传奇（“步骤”，“起点”，“最后一点”）
数字
asolm = sol.a;
NASOLM = SQRT（SUM（ASOLM。^ 2,2））;
情节(nasolm“rx”）
Xlabel（“时间步”）
ylabel（“常态（加速）”）
结束
此代码创建 Fvalt. 函数，它由 findT ．
函数Fval = fvalT (T,选择)
p0 =(-1000、-800、1200);
v0 = (100, -40);
tprob = makeprob（t，p0，v0）;
选择= Optimoptions（选择，“显示”，“关闭”）;
[~, Fval] =解决(tprob选项=选择);
结束
此代码创建 findT 函数。
函数tmin = findt（选择）
DISP（“解决......”）
Tmin = fminbnd (@ (T) fvalT (T,选择),20岁,60);
DISP（“完成”）
结束
版权所有2021 The MathWorks, Inc。
|
类别:
新特性,
优化
  • 打印
  • 发送电子邮件

另请参阅

评论

要发表评论，请点击此处