fmincon with handle function input

Hello, I have a question about the usage of the fmincon routine. I've read the documentation, and the call requires an explicit function with a certain number of variables and coefficients that determine equality and inequality constraints. I wanted to ask if, instead of an explicit function, it's possible to input a handle function, such as a numerical integrator (for example, one that varies initial conditions to minimize the integration output). In the documentation, the only example that shows the use of a handle function still includes an explicit function within the handle.

[Torsten]

Give us an example. It's not clear what you mean.

The functions in which you prescribe the expression you want to minimize and the nonlinear constraints can either be passed to fmincon as function handles (suited for short "one-liners") or functions (suited for difficult computations).

[Hynod]

sorry to be late.

I have the following handle function :

function dx = ODE_1_LanderOrbit(t,x)

    % Assigning constants and computing a_sr perturbation
    r_s   = cspice_bodvrd('Sun', 'RADII', 3);
    r_s   = r_s(2);
    r_m   = cspice_bodvrd('Mercury', 'RADII', 3);
    r_m   = r_m(2);
    mu_m  = cspice_bodvrd('Mercury', 'GM', 3);
    grav  = - mu_m / (norm(x(1:3),2))^3;
    p_sm  = cspice_spkpos('MERCURY', t, 'J2000', 'NONE','SUN');
    p_sl  = x(1:3) + p_sm;
    T_s   = 5800;
    sigma = 5.67e-8 ;
    a_sr  = sigma*T_s^4*( r_s / (norm(p_sl,2)) )^2 * ( p_sl / (norm(p_sl,2)) ) / (3*10^8) * 3 / 3490 / 1000;
    
    % Expressing the differential system
    Grav = [  0, 0, 0, 1, 0, 0;
              0, 0, 0, 0, 1, 0;
              0, 0, 0, 0, 0, 1;
           grav, 0, 0, 0, 0, 0;
           0, grav, 0, 0, 0, 0;
           0, 0, grav, 0, 0, 0];
    dx = Grav*x + [0 0 0 a_sr']';
    
end

This function I've just provided works with ode113 and provides the complete integration of the orbit and consequently the final state, starting from an initial state. I was wondering if it's possible in some way to input the integration itself into fmincon so that fmincon can vary the initial state to minimize, for example, the norm of the final state.

[Sam Chak]

Hi @giacomo

If you want to minimize the norm of the states by varying the initial states, isn't it true that, logically, the initial states beginning from the equilibrium point will give the lowest possible norm value?

tspan  = [0 10];
init   = [0 0];     %  initial states (equilibrium)
[t, x] = ode45(@odefcn, tspan, init);
nx     = norm(x)

nx = 0

% system
function dxdt = odefcn(t, x)
    A    = [0 1; -1 -sqrt(3)];
    dxdt = A*x;
end

[Hynod]

reaching the minimum norm value was an example to let me explain about the issue, there are a lot more things that i have to minimize, so i asked for fmincon routine. Despite of that, i think is not possible to find an equilibrium point because i have no explicit function that computes the final state, it is an integration of a nonhomogeneous system of differential equations. It is true that there are a lot of local minimum of the norm during the orbit, but i want to minimize some values related to the orbit, at the end of the t_span integration time. If it can help, i am talking about variations of orbital parameters to minimize, but my intentions are to understand how fmincon really works with an handle function input of this type.

[Torsten]

p = fmincon(@fun,...)
    
function obj = fun(p)
  [t,x] = ode113(@ODE_1_LanderOrbit,...)
  ...
  obj = ...;
end
  
function dx = ODE_1_LanderOrbit(t,x)
  ...
end

[Hynod]

thanks a lot @Torsten, at least this following script seems to run.

function obj = fun(p)
cspice_furnsh('Kernels_all/mykernelsMERCURY.furnsh')
mu_m = cspice_bodvrd( 'mercury', 'GM', 3) ;
t_i1      = cspice_str2et('2041 FEB 10  02:00:00 UTC');
t_f1      = cspice_str2et('2041  AUG 07 08:01:31 UTC');   
t_span1   = t_i1:30:t_f1;
options = odeset('RelTol', 1.e-16, 'AbsTol', 1.e-16);
[t,S_1] = ode113(@ODE_1_LanderOrbit, t_span1 , p ,options) ;
obj = norm(S_1(1:3));
end

(calling it as you suggested with

p = fmincon(@fun,s_0)

where s_0 is already defined in the workspace )

The problem is that it is not an explicit function, so it seems that fmincon never stops finding the minimum, how can I limitate the iterations? if I am not wrong , I think fmincon calls every time my integrator ODE_1_LanderOrbit searching the minimum value of the integration output varying the initial state, so how to stop it finding the minimum?

[Torsten]

S_1 is a numel(tspan1) x 3 matrix. obj = norm(S_1(1:3)) computes sqrt(S_1(1,1)^2+S_1(2,1)^2+S_1(3,1)^2); . I doubt that this really what you want.

Maybe you mean norm(S_1(end,:)), but I'm not sure.

[Sam Chak]

Hi @giacomo,

Out of curiosity, what does finding a set of initial states, S(0), that minimize the norm(S), physically mean for the lander? Understanding this is crucial in designing the objective function, which requires a blend of art and science—both inherent in human attempts to comprehend and describe the world around us.

From Wikipedia: A lander is a spacecraft that descends towards, then comes to rest on the surface of an astronomical body other than Earth.

Image: Chandrayaan-2 lander.

[Hynod]

@Sam Chak

I am writing a thesis about a spacecraft which lands at the poles of Mercury. The objective is to approach the planet thanks to the perturbation of the sun, to minimize the propellant which is necessary to let the orbit be appropriate to land. So i'am trying to find the initial state which let the orbit approach in the better way the planet. The orbit is not keplerian so it changes every second and finding the initial state which consent to evolve the orbit in the way i want to, is fundamental. About the minimization of the norm, well the spacecraft has to reach a very low distance from the planet, so with some constraints (like the norm has to be greater than the radius of Mercury), i can compute a position of the spacecraft at a very low latitude before the landing approach. (sorry for my english, I am italian :) )

[Hynod]

@Torsten I meant the distance between the center of the planet and the spacecraft at the end of the integration, so it is S_1(end,1:3) , i did't copy that well in the question :) .

[Sam Chak]

*The objective is to approach the planet thanks to the perturbation of the sun, to minimize the propellant which is necessary to let the orbit be appropriate to land.*

Based on the description, this sounds more like an orbit transfer problem than a lander descent problem. The spacecraft needs to be injected into an orbit where you think it is the best to trigger the lander descent event. If I'm not mistaken, the motion of the spacecraft needs to account for the motion of Mercury. However, is designing the objective function as a norm on the states S sufficiently justified for minimizing the burning of propellant?

[Hynod]

@Sam Chak I appreciate that you want to enter more in deep with the problem, so I will be more clear and specific to let you understand that well.

You are rigth, fmincon it is only used to resolve only one of the 4 phases of the space mission : the de-orbiting (you called orbit transfer, but that's ok). The maneuver proposed in the thesis involves the solar perturbation causing the initial orbit (determined by the initial state) to evolve in such a way that the final orbit has a very low pericenter (the closest distance of the spacecraft from the planet in an elliptical orbit), for example, at an altitude of 30 kilometers. Your question is pertinent, and I also considered it during the project phase. The issue is actually related to the fact that if the pericenter were higher (resulting in a relatively high final state norm), it's true that the orbital velocity would be lower, and therefore, the propellant required to slow it down would be less. However, it's also true that during freefall when the spacecraft accelerates radially towards Mercury, it gains additional velocity that will still need to be compensated during landing using additional thrusters. The mission is thus divided into : DE-Orbiting , solid rocket motor stage , coasting\freefall , liquid rocket engine stage.

In a standard space mission, the de-orbiting phase is supported by engines or atmospheric drag , but ''unfortunately'' Mercury has no atmosphere :( .

If you are interested , the reason why i said that there are other parameters that i have to minimize ,is that the spacecraft has to land in a specific area with a specific direction of provenience, so i have also to minimize the angle between the velocity and the prefered direction.

[Walter Roberson]

@Hynod

Mercury has an exosphere . Depending how extensive it is, it might need to be taken into account, https://www.space.com/18644-mercury-atmosphere.html as I suspect it might be enough to provide some heating during landing.

[Hynod]

@Walter Roberson exosphere is a region where the density is approximately zero, the heat you are talking about is refered to the kinetic definition of temperature. When we says that the temperature is hundreds of celcius degrees above the atmosphere , is not the medium temperature in a volume, it is related to the energy of the single particles. So atmospheric drag can't be used and the spacecraft is not heated by the passage in the exosphere. The principal fonts of heat are the radiation form the sun and from mercury itself, which irradiates a lot on the brigth side.

[Walter Roberson]

is the density approximately zero though? It is enough to create bow shocks. https://www.nature.com/articles/s41467-020-18220-2

I have no idea what density would result in significant thermal heating of a landing probe, but I don't think it should be entirely neglected. At the very least I would suggest that measurements of heat and magnetism other factors during the descent could be interesting science.

[Hynod]

The bow shock to which the article refers, is the bowshock generated by the solar wind interacting with the magnetic field of mercury, nothing to do with the bowshocks like those experienced by rockets during their ascent to orbit, and even in that case, the bowshocks occur well before the exosphere . The professor indeed would have let me aware of something so relevant fo the de-orbiting, and MESSENGER orbit follows very fidelty the itegration of my ODE system with the perturbations I consider : 3D body of the Sun , J2 J3 J4 gravitational effects , solar radiation pressure.

[Walter Roberson]

The professor indeed would have let me aware of something so relevant fo the de-orbiting,

Sorry, but No. Professors are not semi-supernatural, able to think of everything. If they were then you would not be writing a thesis: you would be writing a lab report.

If there were trajectory resistence due to shock formation in the exosphere of Mercury, how much would there have to be to affect the trajectory (to beyond the round-off error of the calculations -- more, for example, than the difference between using ode45() and ode78()) ? If there were heating due to the shock, how much would there have to be to be of concern? Now is there evidence that the actual density of the exosphere of Mercury is definitely below those levels, beyond any reasonable margin of error in the current modeling of the exosphere? At the very least your thesis should include a paragraph as to why you can be certain that it was definitely unnecessary to model those effects.

[Hynod]

@Walter Roberson I told you that those shocks you were referring to , were generated by the solar wind interacting with the Mercury's magnetic field. I also told you (believe me or not) that MESSENGER real data from NASA about the orbit, follow very faithfully the model I structured with ODE , so I don't know why you are so insistent. The cspice routine i use are from the S.P.I.C.E toolkit of NASA JPL , is a very attendible font of data about trajectories of spacecraft of the past.

[Yousef]

Yes, you can definitely use a function handle as an input to `fmincon`, and the function that the handle refers to can be quite complex. In fact, many optimization problems involve non-trivial functions, including numerical integrators.

Let's break down your question a bit:

1. **Function Handle for Objective Function**: Your objective function can be a numerical integrator or any other complex routine. Here's a simple way to set it up:

matlab code:

function cost = myObjectiveFunction(x)
% Here, x is the vector of variables you're optimizing.
% You can use x as initial conditions or parameters for your numerical integrator.

% ... Your complex routine, e.g., numerical integration
result = myIntegrator(x); % This is just a hypothetical function

% The objective is to minimize 'result'
cost = result;
end

% Call fmincon
x0 = [initial_guesses]; % replace with your initial guesses
[x_optimal, fval] = fmincon(@myObjectiveFunction, x0, ...);
    ```

2. **Equality and Inequality Constraints**: If you have constraints that involve complex operations or other integrations, you can use function handles for those as well.

matlab code:

function [c, ceq] = myConstraints(x)
% c(x) <= 0 constraints
c = ...; % your inequality constraints
    
% ceq(x) = 0 constraints
ceq = ...; % your equality constraints
    end

% Call fmincon
options = optimoptions('fmincon','SpecifyConstraintGradient',false);
[x_optimal, fval] = fmincon(@myObjectiveFunction, x0, [], [], [], [], [], [], @myConstraints, options);
```

3. **Using Integrators within the Optimization Routine**: If you're using MATLAB's built-in numerical integrators (like `ode45`, etc.), you can easily embed them inside your objective function or constraints:

matlab code:

    function cost = myObjectiveFunction(x)
        [~, Y] = ode45(@(t,y) myODE(t, y, x), tspan, initial_conditions);
        % Process Y to get the cost, if necessary
        cost = ...; % some function of Y
            end
        ```

Here, `myODE` would be another function that defines your differential equations. The point is, `fmincon` doesn't care how the function value is computed. It just needs a function handle to obtain values for given inputs.

In essence, as long as you can define your problem in terms of functions (handles) that `fmincon` can call to get objective values, constraints, and optionally gradients, you can use any complex logic within those functions, including numerical integrators.

[Walter Roberson]

In essence, as long as you can define your problem in terms of functions (handles) that `fmincon` can call to get objective values, constraints, and optionally gradients, you can use any complex logic within those functions, including numerical integrators.

Both fmincon() and ode45() have the restriction that the functions being invoked by them must have continuous first derivatives (and second derivatives too in some cases). Therefore you cannot use any complex logic in the functions.

The particular system being modeled has four phases. If there is a smooth (2 second derivatives) transition between each of the phases, then fmincon() and ode45() should be able to handle it, but if there is an abrupt transition anywhere in the mix then you would need to break up the logic into multiple parts.

Answer 1

[Sam Chak]

HI @giacomo,

I'm unsure if this is what you want to minimize using fmincon(). This is just a simple example without constraints. It returns the solution that is very to the true solution (the equilibrium point).

% search the initial states that minimize the norm using fmincon
initial_guess = [10 -10];
[x0sol, fval] = fmincon(@objfun, initial_guess)

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.

x0sol = 1×2

1.0e-10 * 0.1119 -0.8979

fval = 1.1892e-10

% objective function <-- solve the Lander ODE with ode113 and find the norm
function J = objfun(x0)
    tspan  = [0 10];
    init   = [x0(1) x0(2)];     %  initial states (equilibrium)
    [t, x] = ode45(@odefcn, tspan, init);
    J      = norm(x);
end

% system <-- you already have this one (the Lander ODE)
function dxdt = odefcn(t, x)
    A    = [0 1; -1 -sqrt(3)];
    dxdt = A*x;
end

Answer 2

[Catalytic]

No, the objective function that you give to fmincon can contain anything, but fmincon's algorithms will assume that the input-to-ouput mapping is differentiable and in some cases twice differentiable.

[Hynod]

Can you please take a look at my comment with the script? Maybe my question now is clear.

Answer 3

[Matt J]

The problem is that it is not an explicit function, so it seems that fmincon never stops finding the minimum

It is possible that fmincon is struggling to find the minimum, but that shouldn't be related to the "explicitness" of the function. There certainly is no requirement that fun(p) be implemented in a single-line formula, if that's what worries you.

Note however, that the solution to an ODE can be a discontinuous (and therefore non-differentiable) function of the initial state. That can be a problem, since fmincon assumes fun(p) to be continuous and differentiable. Note also that numerical ODE solvers like ode113() have certain fragilities. This is why MathWorks has posted some relevant guidelines about problems involving ODEs, which you should read.

so it seems that fmincon never stops finding the minimum, how can I limitate the iterations?

You can limit the number of iterations by setting the MaxIterations option,

opts=optimoptions('fmincon','MaxIterations',100);
p=fmincon(fun,____,opts)

Generally, though, if you have to cut fmincon off after some number of iterations, it means the optimization is failing.

[Ajay]

nice

[Hynod]

@Matt J opts was accepted from the routine, but it seems that it integrates more times than i selected in optimoptions , i can say it because at every whole integration the warning that RelTol has been increased appears (this message appears every time i integrate with ode113 and options so there's noting to worry about that) and the message appears more than 3 times (the number of iterations i selected).

[Matt J]

Yes, the number of iterations is not the same as the number of function evaluations. That, you would control with MaxFunctionEvaluations.

But note that because you are not supplying a gradient calculation (SpecifyObjectiveGradient=false), fmincon needs to call your objective fun at least N+1 times to make a finite difference approximation of the gradient. Here, N is the number of unknown variables..

[Hynod]

I don't want to let you waste your time, but it's the first time i use fmincon. If i have this

clc,clear

opts1 = optimoptions('fmincon', 'MaxFunctionEvaluations', 3);
% Compute the values first and store them in double variables
x1 = 4000;
x2 = 0.8;
x3 = ((90 - 1.5) + (90 - 62.09)) * 2 * pi / 360;
x4 = 0;
x5 = ((90 - 5.2) * 2 * pi / 360);
x6 = 0;
x7 = 0;

initial_guess = [x1, x2, x3, x4, x5, x6, x7];

obj = fmincon(@fun, initial_guess, [], [], [], [], [], [], [], opts1);

it continues to integrate more than 3 times, and I don't know why.

Additionally, I don't know why thw initial guess and the final value obj1 are the same (obj is not in the workspace and I don't know why)

function obj1 = fun(p)

cspice_furnsh('Kernels_all/mykernelsMERCURY.furnsh')

mu_m = cspice_bodvrd( 'mercury', 'GM', 3) ;

t_i1      = cspice_str2et('2041 FEB 10  02:00:00 UTC');

t_f1      = cspice_str2et('2041  AUG 07 08:01:31 UTC');   

t_span1   = t_i1:30:t_f1;

options = odeset('RelTol', 1.e-16, 'AbsTol', 1.e-16);
p=cspice_conics([p mu_m]',2);
[t,S_1] = ode113(@ODE_1_LanderOrbit, t_span1 , p ,options) ;
KAPPA=cspice_oscltx(S_1(end,:)',2,mu_m);
  obj1 =norm( [ KAPPA(1:6,1) ; KAPPA(9,1) ; KAPPA(10,1) ] - 10e+5 * [0.023884853422025  0.00000898867626   0.000020223757574  0.000000648481377  0.000018125538251  0.00006273268635   0.000058697621847  0.236174160800096]' ) ;

[Walter Roberson]

fmincon() asked to process 7 variables will immediately call the function 7+1=8 times in order to create a starting point for minimization. Those 7+1 = 8 will be counted in the final count, but the maximum iterations will not be checked until after those 7+1 = 8 calls are done.

[Sam Chak]

Additionally, I don't know why thw initial guess and the final value obj1 are the same (obj is not in the workspace and I don't know why)

Your chosen objective function is the Euclidean norm of the vector KAPPA, κ and it is defined by

KAPPA = cspice_oscltx(S_1(end,:)', 2, mu_m);

If κ directly depends on final state vector and you are trying to find the initial state vector that minimizes , isn't it logically true that fmincon will return ?

Are you trying to determine when to start 'deorbiting' the spacecraft by finding the initial state vector that minimizes the burning of propellant? If so, I think that the objective function needs to be modified to address this as well.

[Hynod]

@Sam Chak I have to minimize the difference between a vector of final orbital parameters chosen by me, and the final vector of orbital parameters computed by the integration , indeed cspice_oscltx routine's output is the vector of the orbital parameters. Unfortunately, fmincon works only with scalar output so I thought that the norm of the vector difference (between the output od cspice_oscltx and the one selected by me) was the best value to minimize, but I am not sure it is :( .

Thanks a lot @Walter Roberson , now I understand why fmincon integrated so many times :)

[Sam Chak]

@Hynod, What do you mean by "fmincon works only with scalar output"?

The objective function needs to be carefully designed. But before that, it is necessary to understand what physical parameter in the system that you want to minimize. For example, if you want to obtain the time-optimal deorbiting maneuver from a non-equilibrium initial state, then the quadratic error function over time can be considered. The Euclidean norm is somewhat related to the quadratic error function. However, you considered only the norm of the final state.

[Hynod]

@Sam Chak I mean, correct me if I am wrong, fmincon can only minimize a scalar value, not an array. I've read that on the documentation so I am quite sure of that but not 100%.

[Walter Roberson]

The objective function for minimization must return a scalar for all minimization algorithms other than gamultiobj and paretosearch

[Hynod]

@Walter Roberson can you tell me the difference between those two routines you have mentioned?

[Sam Chak]

Hi @Hynod,

As you can see in my example below, we can work on a vector (state x) and return a scalar representing the Euclidean norm. This scalar can be used as the objective function value. However, in your code, you take the norm of the final state value, which also returns a scalar. Nevertheless, there is a significant difference in the physical meaning between the vector norm and the end-state norm.

% solving the 1st-order ODE
[t, x] = ode45(@(t, x) - x, [0 10], 1);

% check the size of vector, x(t)
Sx = size(x)

Sx = 1×2

53 1

% vector-based norm
Nx = norm(x)

Nx = 2.3371

% end-state norm
Ne = norm(x(end))

Ne = 4.5600e-05

% exponential decay solution
plot(t, x, 'linewidth', 1.5), grid on

[Walter Roberson]

@Hynod

Unfortunately at the moment I do not really understand the difference between gamultiobj and paretosearch . It looks like they use different algorithms, with gamultiobj() using genetic algorithm approaches, and paretosearch using pattern searching; https://www.mathworks.com/help/gads/paretosearch-algorithm.html

Answer 4

[Hynod]

@Matt J @Catalytic @Torsten @Walter Roberson @Sam Chak I thank you all for the support, fmincon reached the minimum I wanted by minimizing the norm of the difference between the orbital parameters vector I wanted, and the vector it has been actually computed by the integration.

[Hynod]

indeed , i would be very curious about how paretosearch and gamultiobj work, because they not seem to require any initial guess and I don't know why.

[Walter Roberson]

gamultiobj() for one internally generates a number of starting points (number == "population size") based upon the upper and lower bound. You can use the options for gamultiobj to supply a particular starting population.