GPOPS2 是 MATLAB 上的一个软件包,用于求解优化控制问题。下载安装见:GPOPS2学习笔记Ⅰ:介绍,下载,安装
问题结构
我们使用 GPOPS2 的目的都是求解各自的最优控制问题,他们往往长得千奇百怪,因此求解的第一步是将这些问题构造成统一的形式。 GPOPS2 采用了一个P阶段最优控制问题形式,当我们希望它能够求解我们手上的问题时,首先需要将我们的问题形式化为这个标准的问题结构。
P-phase 最优控制问题 的一般形式
定义状态量:,控制量:
,初始时间(or 自变量):
,终端时间(or 自变量):
,积分:
。
代表阶段,
是静态参数。P-phase最优控制问题形式如下:
- 目标函数:
- 动力学约束:
- 事件约束:
- 路径约束:
- 积分约束:
where
GPOPS−II 组织结构
GPOPS−II的组织结构如下,它是用户指导程序如何按照此结构执行,也就是我们具体代码要填充的东西。为了指定要解决的最优控制问题,用户必须编写以下MATLAB函数:(1)an endpoint function;和(2)a continuous function。
endpoint function定义了问题中任何阶段的起点和/或终点、问题中任意阶段的积分和静态参数如何相互关联。端点函数还定义了需要最小化的目标函数。continuous function定义了问题任何阶段动力学的演化,计算问题任何阶段任何积分所需的被积函数,以及问题任何阶段的任何路径约束。接下来,用户必须指定以下变量的下限和上限:
(1) the time at the start and terminus of a phase;
(2) the state at the start of a phase, during a phase, and at the terminus of a phase;
(3) the control during a phase;
(4) the path constraints
(5) the event constraints;
(6) the static parameters.
用GPOPS−II构造一个最优控制问题
下面介绍如何在MATLAB中使用GPOPS−II构造并求解一个最优控制问题。
输入
当我们需要调用 GPOPS2 时,只需在代码中直接调用:
output = gpops2(input)
其中输入是包含关于要解决的最优控制问题的所有信息的用户定义结构体,并且输出是包含通过解决最优控制问题而获得的信息的结构体。
用户定义的结构体 setup(input) 包含必需字段和可选字段。结构体中必需的字段如下:
- name:一个没有空格的字符串,包含问题的名称;
- functions:包含continuous function和endpoint function名称的结构;
- bounds:包含问题中不同变量和约束的下限和上限信息的结构;
- guess:包含对问题中的时间、状态、控制、积分和静态参数的猜测。
可选字段如下:
- auxdata:一种包含辅助数据的结构,可由问题中的不同函数使用。加入auxdata字段消除了在问题中指定全局变量的需要。
- derivatives:指定NLP求解器使用的导数近似值和NLP求解程序使用的导数阶数(“第一阶”或“第二阶”)的结构。
- 网格:指定要使用的网格细化方法类型、网格细化精度公差以及初始网格的信息的结构。
- nlp:指定要使用的nlp解算器以及要在选定的nlp求解器中使用的选项的结构。
- displaylevel:一个整数,取值为0、1或2,并提供在GPOPS−II执行期间发送到MATLAB命令窗口的输出量。
输出
GPOPS−II执行的输出同样是结构体,其中输出包含以下字段:
- result:包含以下字段的结构体:
- solution:每个阶段的最佳时间、状态和控制,以及静态参数向量的最佳值。最优时间、状态和控制分别存储在 solution.phase(p)、time、solution.phase(p).state 和 solution.phase(p).control字段中,而静态参数存储在字段 solution.parameter 中;
- objective:最优控制问题的目标函数的最优值;
- result.setup:GPOPS−II生成结果的 result 的 setup 结构体(the setup structure that produced the result found in result with GPOPS − II;);
- result.nextsetup:如果上一个网格的网格细化过程没有终止(由于未满足网格细化精度公差),则将用于下一个网格细化迭代的设置结构。
- meshhistory:求解NLP的每个网格的解和误差估计(仅当使用网格细化时);
- meshiterations:GPOPS−II进行的网格细化迭代次数(仅当使用网格细化时);
具体每一个字段的详细内容将放在后面的内容中。接下来用一个例子来说明这些字段都是如何填写的。
注意:任何不适用于问题的字段(例如,没有路径约束的问题)都应该完全省略。
例子
BrysonDenham
这个是官方自带的 BrysonDenham 问题的例子,可以看到在调用GPOPS2之前,填写了上述的setup字段作为input。
BrysonDenhamMain.m
%------------------------- Bryson-Denham Problem -------------------------%
% This problem is taken from the following reference: %
% Bryson, A. E. and Ho, Y-C, "Applied Optimal Control: Optimization, %
% Estimation, and Control," Hemisphere Publishing, 1975. %
%-------------------------------------------------------------------------%
clear all
close all
clc
setup.name = 'BrysonDenham';
setup.functions.continuous = @BrysonDenhamContinuous;
setup.functions.endpoint = @BrysonDenhamEvents;
setup.nlp.solver = 'snopt';
setup.derivatives.supplier = 'sparseCD';
setup.derivatives.derivativelevel = 'second';
%setup.scales.method = 'automatic-bounds';
setup.method = 'RPMintegration';
setup.mesh.method = 'hp';
setup.mesh.tolerance = 1e-6; % default 1e-3
N = 1;
setup.mesh.phase.colpoints = 4*ones(1,N);
setup.mesh.phase.fraction = ones(1,N)/N;
%--Problem Setup--
%state 1
x10 = 0;
x1f = 0;
x1min = 0;
x1max = 1/9;
%state 2
x20 = 1;
x2f = -1;
x2min = -10;
x2max = 10;
%state 3
x30 = 0;
x3min = -10;
x3max = 10;
%control
umin = -10;
umax = 5;
%time
t0 = 0;
tmax = 50;
%--Initial Guess--
x1_guess = [x10; x1f];
x2_guess = [x20; x2f];
x3_guess = [x30; x3max];
u_guess = [-10; 0];
t_guess = [t0; tmax];
%--Bounds & Guess--
%--State, Control, Time--
setup.bounds.phase.initialstate.lower = [x10, x20, x30];
setup.bounds.phase.initialstate.upper = [x10, x20, x30];
setup.bounds.phase.state.lower = [x1min, x2min, x3min];
setup.bounds.phase.state.upper = [x1max, x2max, x3max];
setup.bounds.phase.finalstate.lower = [x1f, x2f, x3min];
setup.bounds.phase.finalstate.upper = [x1f, x2f, x3max];
setup.bounds.phase.control.lower = [umin];
setup.bounds.phase.control.upper = [umax];
setup.bounds.phase.initialtime.lower = t0;
setup.bounds.phase.initialtime.upper = t0;
setup.bounds.phase.finaltime.lower = t0;
setup.bounds.phase.finaltime.upper = tmax;
%--No Path Constraints--
%--Guess--
setup.guess.phase(1).state = [x1_guess, x2_guess, x3_guess];
setup.guess.phase(1).control = u_guess;
setup.guess.phase(1).time = t_guess;
%setup.auxdata = 10;
%--Call the Solver--
output = gpops2(setup);
solution = output.result.solution;
figure(1)
plot(solution.phase(1).time,solution.phase(1).state,'-o', 'markersize', 7, 'linewidth', 1.5);
grid on
figure(2)
plot(solution.phase(1).time,solution.phase(1).control,'-o', 'markersize', 7, 'linewidth', 1.5);
grid on
figure(3)
plot(solution.phase(1).time,solution.phase(1).costate,'-o', 'markersize', 7, 'linewidth', 1.5);
grid on
BrysonDenhamContinuous.m
function phaseout = BrysonDenhamContinuous(input)
%t = input.phase.time;
%x1 = input.phase.state(:,1);
x2 = input.phase.state(:,2);
%x3 = input.phase.state(:,3);
u = input.phase.control;
dx1 = x2;
dx2 = u;
dx3 = (u.^2)./2;
phaseout.dynamics = [dx1, dx2, dx3];
BrysonDenhamEvents.m
function output = BrysonDenhamEvents(input)
x3f = input.phase.finalstate(3);
%--Cost--
output.objective = x3f;