# Using fmincon to solve objective function in integral form.

22 views (last 30 days)

Show older comments

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))

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 Comments

Dyuman Joshi
on 10 Apr 2024

Bruno Luong
on 10 Apr 2024

### 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 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)

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!