AI与控制（二）从优化到最优控制，从动态规划到强化学习--1

优化问题，尤其静态优化问题，在控制系统设计中随处可见，例如基于燃油经济性和驾驶体验的多目标优化的汽车发动机 MAP 标定，基于性能指标优化的飞行器结构设计参数优化，以实验数据与模型输出匹配为目标的电池 RC 等效电路模型标定等等，他们都是通过构建目标函数（某个值的性能最大，或者某两个值之间的差距最小），然后调用优化算法实现设计变量寻优。

当设计变量变成一个函数，而这个函数对系统指标的影响又在时间先后上取决于一个动力学约束，这时候我们可能依然可以通过离散的方式将问题变成静态优化问题，当然这类问题也可以通过最优控制理论来实现，利用庞特里亚金极小值原理，动态规划来求解。例如混动车辆的 ECMS 算法，输出电池和发动机的能量分配序列满足电池 SOC 平衡的基础上油耗最低。

随着 AI 的引入，即使我们对于系统或模型一无所知，我们又可以通过试错的方式来获得一个长期奖励较优的控制器，用于处理序列决策问题，例如自动驾驶车辆或智能机器人的控制器控制序列，也就是强化学习的思路。本文接下来通过MATLAB示例来简单介绍这些概念的思想。

优化问题

对于一个普通的静态优化问题，可以描述为求解最优变量 x 使得 f(x) 最小[1]，

同时满足约束条件：

其中 x 是 n 维的设计变量向量，f(x) 是目标函数（通常是标量），和 m 维的函数向量 G(x),用于计算 x 处的等式和非等式约束。例如，求解这个示例：

目标函数

约束

这个示例是典型的包含非线性约束的平滑问题，为了迭代高效，可以在目标函数和约束函数计算时中同时给出梯度。

function [f,gf] = onehump(x)

% ONEHUMP Helper function for Tutorial for the Optimization Toolbox demo

% Copyright 2008-2009 The MathWorks, Inc.

r = x(1)^2 + x(2)^2;

s = exp(-r);

f = x(1)*s+r/20;

gf = [(1-2*x(1)^2)*s+x(1)/10;

-2*x(1)*x(2)*s+x(2)/10];

end

function [c,ceq,gc,gceq] = tiltellipse(x)

% TILTELLIPSE Helper function for Tutorial for the Optimization Toolbox demo

% Copyright 2008-2009 The MathWorks, Inc.

c = x(1)*x(2)/2 + (x(1)+2)^2 + (x(2)-2)^2/2 - 2;

ceq = [];

gc = [x(2)/2+2*(x(1)+2);

x(1)/2+x(2)-2];

gceq = [];

end

构建完目标函数和约束函数，这样就可以调用 fmincon 求解器求解这个非线性平滑约束的优化问题，

options = optimoptions(options,...

'SpecifyObjectiveGradient',true,...

'SpecifyConstraintGradient',true);

[x,fval,exitflag,output] = fmincon(@onehump,x0,[],[],[],[],[],[], ...

@tiltellipse,options);

xold = x

xold = 2x1

- 0.9727

0.4686

于是得到约束条件下的最优解（红色的点）:

这类有约束的问题中，我们的优化变量是数值向量，我们可以通过 MATLAB 内置的求解器来求解极值点（算法参考[2]），也可以使用符号数学工具箱，利用函数的极值条件，手动构建 Lagrange multiplier（将等式约束问题转换为无约束问题）或者KKT条件进行求解，例如，针对 Lagrange multiplier 方法，在 g(x)=0 的等式约束条件下求解 f(x) 最小，可以先构造 Lagrangian 函数：

需要满足极值条件：

我们可以用这种方法求解优化问题（1-1）（因为有不等式约束，下面的脚本添加了 s1 项变成等式约束，读者也可以尝试使用 KKT 条件进行含不等式约束问题的求解），可以得到和数值求解相同的结果。

syms x y s1 lambda

f = x.*exp(-x.^2-y.^2)+(x.^2+y.^2)/20;

g = x.*y/2+(x+2).^2+(y-2).^2/2-2+s1^2==0;

L = f + lambda*lhs(g);

dL_dx = diff(L,x)==0;

dL_dy = diff(L,y)==0;

% derivative of L with respect to lambda

dL_dlambda = diff(L,lambda) == 0;

dL_ds1=diff(L,s1) == 0;

% build the system of equations

system = [dL_dx; dL_dy; dL_dlambda;dL_ds1];

% solve the system of equations

[x_val, y_val, lambda_val, s1_val] = ...

solve(system, [x y lambda s1], 'Real', true);

% show results in a vector of data type double

results_numeric = double([x_val; y_val])

results_numeric = 2x1

- 0.9727

0.4686

变分法求泛函极值

现在如果我们的待优化变量不是数值向量，而是一个函数，例如这样一个最简变分问题，J(x) 是一个泛函，它是函数 x(t) 的函数，我们希望求解最优的函数 x*(t), 使得性能指标 J(x) 取得极小值。

其中 t1, t2 是常数，x(t) 是连续二阶可微函数, L(t,x(t),Ẋ(t)) 对其所有变量连续二阶可微。这个问题就是一个求解泛函极值的问题。类比于静态优化求解 f(x) 极值条件时，通过给 x 一个扰动 x+△x 来计算一阶导数 ▽f/▽x 让其等于零从而获得极值条件，利用相似的思路，对于自变量 x(t) 是函数的情况我们此处取任意一个可微扰动函数 η(t) 在边界 t1 和 t2 处的值为零，如果 x(t) 是极值点，有 J(x(t))≤J(x(t)+ε*η(t)) ，其中 ε*η(t) 就叫函数 x(t) 的变分，记为 δx(t)，变分引起的泛函 J(x) 的变化即泛函增量：△J(x,δx)=J(x(t)+δx(t))-J(x(t)), 类比于函数增量的泰勒展开 △y=△xT▽y(x)+o(||△x||) 的函数增量可以写成 “变分” △x 的线性函数与高阶无穷小加和的形式，泛函增量 △J(x,δx) 也可以写成变分 δx 的线性泛函 δJ(x(t),δx) 与高阶无穷小加和的形式 △J(x,δx)=δJ(x(t),δx)+o(||δx||) 。J(x(t)+ε*η(t)) 可以看作 ε 的函数,记 ψ(ε)=J(x(t)+ε*η(t)), ψ(0)=J(x(t))[4] 于是有：

其中

是泛函 J(x) 对函数 x(t) 的泛函微分[3]。如果泛函 J(x) 取得极值，则有 δJ(x,δx)=0，通过展开推导以及 η(t) 的任意性，可以得到极值的必要条件：欧拉-拉格朗日方程，即

我们使用一个示例演示利用 MATLAB 求解泛函微分，从而求得极值函数。质量为 m 的质量块与弹簧常数为 k 的弹簧相连，求它的欧拉-拉格朗日方程。定义系统的动能 T，势能V，和拉格朗日量 L。拉格朗日量是动能和势能之差。

syms m k x(t)

T = 1/2*m*diff(x,t)^2;

V = 1/2*k*x^2;

L = T - V

L(t) =

在拉格朗日力学中，系统的作用泛函等于拉格朗日函数随时间的积分，

欧拉-拉格朗日方程描述了 S[x(t)] 是平稳的系统的运动。所以求被积函数L的泛函导数并令其等于0。

eqn = functionalDerivative(L,x) == 0

eqn(t) =

eqn 是描述质量-弹簧振荡的微分方程。用dsolve求解eqn。假设质量m和弹簧常数k为正。设振荡幅度的初始条件为 x(0)=10; 初始速度为0

assume(m,'positive')

assume(k,'positive')

Dx(t) = diff(x(t),t);

xSol = dsolve(eqn,[x(0) == 10, Dx(0) == 0])

xSol =

当然，相较于最简变分问题（2-1），实际中，我们要求解的系统通常是有约束的，例如，系统有一个动力学方程：

利用 Lagrange multiplier 方法，可以构造增广泛函

从而可以继续利用 Euler–Lagrange 方程得到含动力学约束问题取得性能极值的必要条件。

庞特里亚金极小值原理 (PMP)

基于变分法，前苏联科学家庞特里亚金证明了解决更一般最优控制问题的必要条件[4]。

最优控制问题：

对于一个动态系统，状态变量 x∈ℝn, 控制变量 u∈u, 已知系统的动力学方程

终端时刻自由或固定，终端状态自由或固定。现在需要求解最优控制函数（轨迹）u:[t0,tf]→u,使得性能指标泛函 J 取得极小值。其中 J 定义为：

其中为终端时刻性能指标，L 为瞬时指标。

PMP 证明了性能指标取得全局极小值的必要条件，系统的动力学约束方程可以通过引入时变 Lagrange multiplier 向量 λ(t),也叫协态变量, 构造生成定义在 t∈[t0,tf] 上的 Hamiltonian 函数：

极值条件：

最优控制 u*(t), 最优状态轨迹x*(t) 以及对应的协态轨迹 λ*(t),对于任意的容许控制 u(t)∈u,在任意时刻 t∈[t0,tf], 都满足

协态方程:

边界条件：

我们参考 FileExchange 资源[5]示例演示利用 PMP 求解二阶系统的终端时刻自由的问题。资源[5]引用了一个二阶系统的最优控制问题的求解。

定义性能指标

给定初始条件和终端条件

我们尝试使用 PMP 原理进行求解。

极值条件

首先构建 Hamiltonian 函数和求解极值条件得到 u(t) 的表达。

syms x1(t) x2(t) p1(t) p2(t) u tf t

% 定义动力学约束方程右侧

DFx1 = x2;

DFx2 = u;

% 计算 Hamiltonian 函数,并根据极值条件获得 u(t) 的表达

H = 0.5*u^2+p1*DFx1 + p2*DFx2;

sol_u = solve(diff(H,u)==0,u);

% 将最优 u(t) 的表达带入 DFx2

DFx2 = subs(DFx2,u,sol_u);

状态方程和协态方程

这里直接可以套用 PMP 的协态方程

% 构建 PMP 的状态方程

deqn1 = diff(x1,t)==DFx1;

deqn2 = diff(x2,t)==DFx2;

% 构建 PMP 的协态方程

deqn3 = -diff(p1,t)==diff(H,x1);

deqn4 = -diff(p2,t)==diff(H,x2);

边界条件

除了状态约束，我们通过 PMP 的边界条件可以知道，对于自由终端状态问题，还需满足

本问题因为 h≡0，而且仅是 x2(tf) 自由，因此有边界条件

% 边界条件

cond1 = x1(0)==1;

cond2 = x2(0)==2;

cond3 = x1(tf)==3;

cond4 = p2(tf)==0;

% 通过边界条件，求得状态/协态和自由时间 tf 的表达

sol = dsolve([deqn1 deqn2 deqn3 deqn4],[cond1 cond2 cond3 cond4]);

% 利用上述求得的最优表达和边界约束条件 H(x(tf),u(tf),λ(tf),tf)=0 计算得到 tf

H_new = subs(H,u,sol_u);

sol_tf = solve(subs(H_new,[p1 p2 x2 t],[sol.p1 sol.p2 sol.x2 tf])==0);

% 其中有一个解是 u=0, 我们去掉这个。

sol_tf = sol_tf(2)

x1 = subs(sol.x1,'tf',sol_tf)

x2 = subs(sol.x2,'tf',sol_tf)

p1 = subs(sol.p1,'tf',sol_tf)

p2 = subs(sol.p2,'tf',sol_tf)

u = subs(sol_u)

除了使用符号计算进行解析求解，还可以参考资源[5]将问题作为二点边值问题（two point boundary value problem）,从而使用 bvp4c[6] 进行数值求解。文中介绍了自由时间和自由终端状态问题的处理方法。

免费分享一些我整理的人工智能学习资料给大家，整理了很久，非常全面。包括一些人工智能基础入门视频+AI常用框架实战视频、图像识别、OpenCV、NLP、YOLO、机器学习、pytorch、计算机视觉、深度学习与神经网络等视频、课件源码、国内外知名精华资源、AI热门论文等。

下面是部分截图，加我免费领取