# Using fmincon to solve objective function in integral form.

12 views (last 30 days)
Joseph Criniti on 10 Apr 2024
Commented: Joseph Criniti on 19 Apr 2024
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 CommentsShow NoneHide None
Dyuman Joshi on 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 on 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.

Sign in to comment.

### Accepted Answer

Bruno Luong on 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 CommentsShow 13 older commentsHide 13 older comments
Bruno Luong on 17 Apr 2024
Edited: Bruno Luong on 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 on 19 Apr 2024
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.

Sign in to comment.

### Categories

Find more on Numerical Integration and Differentiation in Help Center and File Exchange

R2020b

### Community Treasure Hunt

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

Start Hunting!