计算科学 - 自动生成轨迹优化约束 - 吾爱随笔录

我有兴趣从上面链接中提到的问题执行轨迹优化。

我想将以下内容作为动态约束提供给 MATLAB 的fmincon.

\frac{d ϕ}{d t} = - M^{T} D M ϕ + W (1)

$\frac{d \phi}{dt} = -M^TDM\phi + W \hspace{1cm}(1)$

$\phi$ 是一个带有变量 [ABCDEF] 的向量。

我一直在参考这篇文章，它说明了梯形搭配是如何完成的

例如，如果是一个标量 $\phi$

\dot{ϕ} = f (2)

$\dot\phi = f \hspace{1cm} (2)$

到对上述内容进行积分，并使用梯形求积进行近似，得到， $t_k$ $t_{k+1}$

ϕ_{k + 1} = ϕ_{k} + \frac{1}{2} h_{k} (f_{k + 1} + f_{k}) (3)

$\phi_{k+1} = \phi_k + \frac{1}{2}h_k(f_{k+1} + f_k) \hspace{1cm} (3)$

这将转化为 k = 0 到 N-1的一组非线性等式约束 ( ref )。

我想知道如何自动生成这些约束。是否可以自动生成这些约束并将其作为句柄传递给fmincon. 因为可以手动编写单个变量，但对于多变量，它变得复杂。

我还查看了可以在 MATLAB 中使用的存储库。但是，我不确定这是否可以用来生成约束。

另外，我想了解方程（3）实际上应该如何写为等式约束。是否应该求解 (1) 中的 ode 系统并使用值来评估等式 (3) 的 RHS？RHS 会发生什么？ $\phi$

对上述问题的任何建议/解释都将非常有用。

编辑：我尝试按照此处给出的示例代码来了解如何实现 ceq（平等约束）

例如，在上面链接中提供的代码中

function [ c, ceq ] = double_integrator_constraints( x )

给出。

有人可以解释输入参数 x 是如何计算的吗？另外，我不明白如何ydesiredend在平等约束中提及。

在我的问题中，只为动态约束定义了初始条件，而终端条件不可用。在那种情况下，我不确定如何定义等式约束。

EDIT2：更新：我可以使用 fminunc 和 lsqnonlin 解决这个问题。

Dhat0 = %input vector 
% fun   = @objfun;
% [Dhat,fval] = fminunc(fun, Dhat0)

%% lsqnonlin
Dhat = lsqnonlin(@(Dhat) objfun(Dhat),Dhat0)


function f = objfun(Dhat)

%% Integrator settings
tspan = %tspan 
options = odeset('abstol', 1e-10, 'reltol', 1e-9);

%% generate exact solution
    phi0 = % initial condition vector
    [t, phi]  = ode15s(@(t,phi) exact(t,phi), tspan , phi0 ,options);


%% generate approximate solution

    [t, phi_tilde]  = ode15s(@(t,phi_tilde) approx(t,phi_tilde, Dhat), tspan , phi0 ,options);


%% objective function for fminunc
    % diff = (phi - phi_tilde).*(phi - phi_tilde);
    % f = sum(diff, 'all')

%% objective function for lsqnonlin
    f  = phi - phi_tilde
end

我仍然有兴趣了解如何解决约束优化问题。我正在尝试以这种方式解决问题。

我想了解noncol应该如何从 fmincon 调用。

nonlcon = @defects;
Dhat= fmincon(@objfun,Dhat0,A,b,Aeq,beq,lb,ub,nonlcon)

对于我的系统，A,b,Aeq,beq,lb,ub = []

但是，我不确定从哪里传递这些论点 defects(dt,x,f)。

function [c ceq] = defects(dt,x,f)
    % ref: https://github.com/MatthewPeterKelly/OptimTraj
    % This function computes the defects that are used to enforce the
    % continuous dynamics of the system along the trajectory.
    %
    % INPUTS:
    %   dt = time step (scalar)
    %   x = [nState, nTime] = state at each grid-point along the trajectory
    %   f = [nState, nTime] = dynamics of the state along the trajectory
    %
    % OUTPUTS:
    %   defects = [nState, nTime-1] = error in dynamics along the trajectory
    %   defectsGrad = [nState, nTime-1, nDecVars] = gradient of defects


nTime = size(x,2);

idxLow = 1:(nTime-1);
idxUpp = 2:nTime;

xLow = x(:,idxLow);
xUpp = x(:,idxUpp);

fLow = f(:,idxLow);
fUpp = f(:,idxUpp);

% This is the key line:  (Trapazoid Rule)
defects = xUpp-xLow - 0.5*dt*(fLow+fUpp);
ceq = reshape(defects,numel(defects),1);
c = [];
end
end

任何建议都会非常有帮助