How to use InfiniteOpt.jl to solve an optimal control with multiple stages?

[ronisbr]

Hi!

I am planning to use InfiniteOpt.jl to find the optimal control for a launcher including its reentry. This problem has many stages and I need to change the dynamics and the constraints at each stage. In the past, I was using PSOPT (C++). I am wondering if this kind of problem can be solved using InfiniteOpt.jl.

[pulsipher]

It should be possible to solve this problem class with InfiniteOpt.jl. If you can provide a basic formulation to exemplify what you have in mind, I can help you determine how to code it up.

[ronisbr]

Hi @pulsipher !

Thanks! I am trying to figure out a simple model that is representative. The main ideia is that we must have phases. Let's suppose we have a rocket. Hence:

• Phase 1: From t = 0 until the height is 200m
• Phase 2: From phase 1 until the fuel is over

If we are in Phase 1, the thrust angle is computed to maximize the distance from the platform without violating some constraints. If we are in Phase 2, the thrust angle is computed to minimize other index.

I know we can use InfiniteOpt.ifelse to have branches in the dynamics. However, I could not find how can I store a integer that indicates the phase.

[ronisbr]

I think the biggest problem here is to change the phase. There is not a variable that I can use to verify in which phase we are. I need some kind of event tracking to switch between the models.

[pulsipher]

The key would be deriving a single NLP formulation where the variables are indexed over a continuous domain. This will likely involve some reformulation tricks in like manner to the minimum time problem in #219. It seems there are essentially two coupled optimization problems being solved here. One naïve way to implement this would be to use a binary variable to indicate which phase (problem) is active at different domain indices and then use big-M constraints; however, this would result in an MINLP. Another way would be to reformulate the problem in terms of a progress parameter to index the variables and then use DomainRestrictions to enforce the proper constraints in each phase based on the height criterium. Unfortunately, it is difficult for me to give more precise advice without more information.

[ronisbr]

Hi @pulsipher !

Thanks for the answer. Looking at PSOPT, I think the correct way here is to create many constraints given the number of phases. Hence, depending on the phase, I set the constraint to 0 or to the correct value. At least, this is my understanding about how PSOPT works when dealing with multiple stages. In my particular case, the state will be changed given the simulation time (I could approximate the problem using this consideration). Hence, it is easy to track what is the current phase.

[pulsipher]

Ok, based on my reading of PSOPT, the multiple phases appear to amount to path constraints enforced over sub-intervals of of the overall time interval. If this understanding is correct, then we can simply enforce the desired constraints in InfiniteOpt.jl using DomainRestrictions as shown below:

using InfiniteOpt

# Define the time horizon info
t1 = 5 # end of phase 1
tf = 10 # final time

# Create the model
model = InfiniteModel()
@infinite_parameter(model, t in [0, tf], num_supports = 100) # an arbitrary discretization choice
@variable(model, 0 <= y <= 1, Infinite(t))

# Add the phase 1 constraint
@constraint(model, phase1, sin(y) <= 0, DomainRestrictions(t => [0, t1]))

# Hack to determine the next discretization point after t1
ts = supports(t)
t_next = ts[findfirst(isequal(t1), ts) + 1]

# Add the phase 2 constraint
@constraint(model, phase2, cos(y) >= 0, DomainRestrictions(t => [t_next, tf]))

# Add a linking constraint between phases if needed
@constraint(model, linking, y(t1) == y(t_next))

Note that the little hack is needed since DomainRestrictions are currently limited to inclusive interval bounds, a rewrite is underway to remove this restriction. I hope this helps.

[ronisbr]

Hi @pulsipher !

Awesome! Thank you very much. This example will help me a lot!

[ronisbr]

Hi @pulsipher !

I have another question. My objective function is an integral (something like integral((ua / ua_max)^2 + (ub / ub_max)^2, t). Where ua and ub are variables. However, the constants ua_max and ub_maxdepend on the current phase. What is the suggested method to take this into account?

[pulsipher]

You can use a parameter_function to define time-dependent constants. For instance:

using InfiniteOpt

# Define the time horizon info
t1 = 5 # end of phase 1
tf = 10 # final time

# Create the model
model = InfiniteModel()
@infinite_parameter(model, t in [0, tf], num_supports = 100) # an arbitrary discretization choice
@variable(model, ua, Infinite(t))
@variable(model, ub, Infinite(t))

# Define ua_max(t) via the do-block syntax
ua_max = parameter_function(t, name = "ua_max") do t
                    if t < t1
                        return 2.0
                    else 
                        return 10.2
                    end
                 end

# Define ub_max(t) using the macro syntax
ub_max_value(t) = t < t1 ? 3.0 : 42.0
@parameter_function(model,  ub_max == ub_max_value(t))

# Define the objective
@objective(model, Min, integral((ua / ua_max)^2 + (ub / ub_max)^2, t))

[ronisbr]

Hi @pulsipher !

First of all, thank you very much for all the help. It has been very simple using InfiniteOpt.jl to solve my problem and everything seems to be working very fine! Congratulations for this amazing package.

I am now facing a new challenge from what I described. In my problem, we have phases (7) as I explained. For the initial phases, where the transition happens in pre-defined time-steps, everything is working perfectly using DomainRestrictions. However, I am facing problems for a phase that starts when a variable starts to decrease (when we reach the maximum altitude). In this case, the time is not defined, but I need to change the dynamics of a variable when this happens.

Do you have any idea what is the best approach here?

EDIT: Just to complement what I am planning as my current approach:

1. Let the infinite parameter t by a continuous variable from 0 to 6. Hence, phase 1 lasts from 0 to 1, phase 2, 1 to 2, and so on.
2. Each phase has a time scale parameter that defines the phase duration.
3. Add the constraints required to obtain the time scales in the phases that we do not know when the change occurs.

The problem here is to define the constraint / objective in those phases that the time is free. For example, in a specific phase, I need to find the time scale that leads to the maximum height. However, I cannot add this to the objective function because it just does not make sense in any other phase.

[pulsipher]

The quick answer would be that for phases that cannot be defined by prescribed times, some reformulation (change of variables is needed) such that the beginning and end of the infinite parameter is known for a particular phase. For instance, in minimum time problems it is possible to scale the problem in terms of the final time (see #188) and in fastest path problems it is possible to reformulate in terms of progress along a path (see #219).

In this case, it might be most convenient to have a different infinite parameter for each phase (using time or some other derived quantity for each) and add the constraints for each phase. The constraints can be linked together using appropriate point variables at the phase boundaries. Note this approach would even remove the need to use DomainRestrictions (and the hack I previously suggested).

[ronisbr]

Hi @pulsipher !

Thanks for the help. In the case I create an infinite parameter for each phase (let's say I have then multiple times), then I need to define variables for each phase? Something like:

@infinite_parameter(model, tp1 ∈ [0, 2])
@infinite_parameter(model, tp2 ∈ [2, 4])

@variable(model, hp1 ≥ 0, Infinite(tp1))
@variable(model, hp2 ≥ 0, Infinite(tp2))

Or is there a better approach?

[pulsipher]

Yes, exactly. In the end, mutli phase problems essentially correspond to multiple coupled optimization problems (coupled via the boundaries and the objective).

[ronisbr]

I will try that, thank you one more time @pulsipher ! :)

[ronisbr]

Hi @pulsipher !

I just want to let you know that this suggestion to have one infinite parameter per phase worked beautifully! I did not finish everything yet, but it really looks promising. Thanks for all the help and for the amazing package!

[ronisbr]

Hi @pulsipher

I just want to let you know that I managed today to implement an optimizer for a very complex, 8 phase problem. InfiniteOpt was amazing and thank you very much for all the advices! The way we write the constraints and expressions in InfiniteOpt makes everything so much clean and easy to understand!

The only thing I miss in Julia ecosystem compared with PSOPT or GSOPT-II is some mechanism to perform mesh refinement for me.

[pulsipher]

Thank you for the update and the kind words! Mesh refinement is on the future wish list for InfiniteOpt, just gotta find the time and right person to do it :)