GPOPS列车最优运控

GPOPS列车最优运控

动力学方程

$$\begin{gathered} \frac{dt}{ds}=\frac{1}{v} \\ \frac{dv}{ds}=\frac{u-r(v)-g(s)}{v} \end{gathered}$$
目标函数
$$J={\int }_{s_0}^{s_f}max(u,0)ds$$
代码，没考虑坡度和限速，所以只有一段

clear clc setup limits guess
tic;
global constants
 
constants.D_total = ;	% 总路程
constants.T_total = ;	% 总时间
constants.A = ;			% 戴维斯系数
constants.B = ;			% 戴维斯系数
constants.C = ;			% 戴维斯系数
constants.M = ;			% 列车质量
constants.PowerMax = ;	% 功率限制
constants.velMax = ;	% 速度限制
constants.uMin = ;		% 控制限制
constants.uMax = ;
 
s0 = 0; sf = constants.D_total;
t0 = 0; tf = constants.T_total;
v0 = 1; vf = 1;
 
% start phase
s0Min = s0; s0Max = s0;
t0Min = t0; t0Max = t0;
v0Min = v0; v0Max = v0;
 
% during phase
tMin = 0; tMax = tf;
vMin = 1; vMax = constants.velMax;
 
% terminus phase
sfMin = sf; sfMax = sf;
tfMin = tf; tfMax = tf;
vfMin = vf; vfMax = vf;
 
uMin = constants.uMin;
uMax = constants.uMax;
 
iphase = 1;
limits(iphase).meshPoints       = [-1 +1];
limits(iphase).nodesPerInterval = 99;
limits(iphase).time.min         = [s0Min  sfMin];
limits(iphase).time.max         = [s0Max  sfMax];
limits(iphase).state.min(1,:)   = [t0Min tMin tfMin];
limits(iphase).state.max(1,:)   = [t0Max tMax tfMax];
limits(iphase).state.min(2,:)   = [v0Min vMin vfMin];
limits(iphase).state.max(2,:)   = [v0Max vMax vfMax];
limits(iphase).control.min      = -uMin;
limits(iphase).control.max      = uMax;
limits(iphase).parameter.min    = [];
limits(iphase).parameter.max    = [];
limits(iphase).path.min         = [];
limits(iphase).path.max         = [];
limits(iphase).event.min        = [];
limits(iphase).event.max        = [];
 
%% guess
guess(iphase).time          = [s0; sf];     % increase
guess(iphase).state(:,1)    = [t0; t0];
guess(iphase).state(:,2)    = [v0; vf];
guess(iphase).control       = [0; 0];
guess(iphase).parameter     = [];
 
 
%% setup
setup.name              = 'Train-Problem';
setup.funcs.cost        = @TrainCost;
setup.funcs.dae         = @TrainDae;
setup.linkages          = [];
setup.limits            = limits;
setup.guess             = guess;
setup.derivatives       = 'finite-difference';
setup.checkDerivatives  = 0;
setup.autoscale         = 'on';
setup.mesh.tolerance    = 1e-1;
setup.mesh.iteration    = 1;
setup.mesh.nodesPerInterval.min = 5;
setup.mesh.nodesPerInterval.max = 10;
 
[output,gpopsHistory] = gpops(setup);
toc;
solution = output.solution;

dae函数

function [dae] = TrainDae(sol)
    global constants
    s = sol.time;
    t = sol.state(:, 1);
    v = sol.state(:, 2);
    u = sol.control;
    tdot = 1 ./ v;
    r = (constants.A + constants.B * v + constants.C * v .^ 2) / constants.M;
    vdot = (u - r) ./ v;
    dae = [tdot vdot];
end

cost函数

function [Mayer,Lagrange] = TrainCost(sol)
    s0 = sol.initial.time;
    u  = sol.control;
    Mayer = zeros(size(s0));
    Lagrange = max(u, 0);
end

控制曲线

速度曲线