Using fmincon to solve objective function in integral form.
15 次查看(过去 30 天)
显示 更早的评论
Hello, I'm attempting to solve an optimal control problem using the fmincon function however I have very little experience with it. I'm wondering if there is some fundamental flaw to my code or if there is a more obscure mistake like using the trapz function for integration. My current batch of code is no longer spitting out errors so trying to troubleshoot my way to a solution has slowed down.
The problem is as follows:
And here is the code:
%% Optimization
clear all
close all
clc
%Initialize constants
global U tspan N div
tspan = 10;
div = 10;
N = tspan*div;
control_ic = [linspace(0, 10, 1/5*N), linspace(10, 0, 4/5*N)]; %Control profile guess
%% Routine
u_t = fmincon(@fun, control_ic, [], [],[],[], [], [], [], optimoptions('fmincon', 'MaxFunctionEvaluations', 1e5))
function [J] = fun(u) %Objective function
global U tspan div
U = u;
dynamic_ic = [2;2]; %Initial Conditions
[t, X] = ode45(@dynamics, [0 tspan], dynamic_ic); %Find x, xd using u
x = X(:, 1);
xd = X(:, 2);
x_cost = trapz(tspan/length(t), 2*x.^2 + 2*xd.^2); %Calculate objective function
finalx_cost = x(length(x))^2;
control_cost = trapz(div^(-1), .1*u.^2);
J = x_cost + finalx_cost + control_cost;
end
%System Dynamics
function [X] = dynamics(t, x)
global U tspan N
xd = x(2);
vd = -.2*x(2) - 3*x(1) + interp1(linspace(0, tspan, N), U, t);
X = [xd; vd];
end
2 个评论
Dyuman Joshi
2024-4-10
Programmatically, I would suggest you to remove the 1st three lines of code, and, remove the global variables and provide them as inputs to the functions and call accordingly.
Bruno Luong
2024-4-10
Also on the control point of view you should regularize the control u by adding a penaty on 2-norm of du/dt in the cost function. This avoid oscillation.
采纳的回答
Bruno Luong
2024-4-10
trapz(tspan/length(t), 2*x.^2 + 2*xd.^2)
This looks wrong to me, you seem to assume t is equi distance (even in that case you should duvide by length(t)-1).
I woild say this formula
trapz(t, 2*x.^2 + 2*xd.^2)
returns what you want. Warning I do not test, and possibly there is still other mistake in your code.
15 个评论
Bruno Luong
2024-4-17
编辑:Bruno Luong
2024-4-17
I extend my linear ode equation solver with the forcing term u(t) - the control - to piecewise linear (it was piecewise constant).
It must be out there someone who wrote formula recipe for linear ode with generic polynomial forcing term. I do it by hand so it is a little bit painful to increase the degree.
EDIT I think have derived a general closed formula to solve linear ode with constant coefficient and general polynomial forcing U(t).
dX/dt = A*X + U(t)
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!