# Using fmincon to solve objective function in integral form.

32 visualizzazioni (ultimi 30 giorni)
Joseph Criniti il 10 Apr 2024
Commentato: Joseph Criniti il 19 Apr 2024 alle 20:06
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))
Local minimum possible. Constraints satisfied. fmincon stopped because the size of the current step is less than the value of the step size tolerance and constraints are satisfied to within the value of the constraint tolerance.
u_t = 1x100
-4.7060 17.2991 2.8918 -25.1503 23.6414 5.3004 -0.3443 -1.4582 8.0192 13.3754 5.3292 -2.3025 11.8676 6.8500 9.7229 6.1461 6.4091 8.8996 12.3205 10.0425 3.6941 14.0695 2.2968 11.0991 9.4470 13.5367 -0.5043 12.9030 5.6669 8.4720
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
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 CommentiMostra NessunoNascondi Nessuno
Dyuman Joshi il 10 Apr 2024
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 il 10 Apr 2024
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.

Accedi per commentare.

### Risposta accettata

Bruno Luong il 10 Apr 2024
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 CommentiMostra 13 commenti meno recentiNascondi 13 commenti meno recenti
Bruno Luong il 17 Apr 2024
Modificato: Bruno Luong il 17 Apr 2024
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)
Joseph Criniti il 19 Apr 2024 alle 20:06
Thank you for your continued support. I appreciate you going out of your way to optimize the code. Your insight has been very helpful for my implementation fmincon in optimal control problems in the future. Once again, I can't thank you enough.

Accedi per commentare.

### Categorie

Scopri di più su Numerical Integration and Differentiation in Help Center e File Exchange

R2020b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by