control-toolbox / CTBase.jl

Fundamentals of the control-toolbox ecosystem
MIT License
7 stars 2 forks source link

Suppress Index in constraints #155

Open jbcaillau opened 2 weeks ago

jbcaillau commented 2 weeks ago

... but keep for (variable) times IMG_3235

PierreMartinon commented 2 weeks ago

Can you elaborate please :D ?

PierreMartinon commented 1 week ago

Hi @ocots @jbcaillau What is the new syntax for a basic initial condition x(t0) = x0 ? The version below is not accepted anymore. I updated the time! call but got stuck on the constraint! (did not find a similar example in the doc)

ocp = Model(variable=true)
state!(ocp, 3)
control!(ocp, 1)
variable!(ocp, 1)
time!(ocp, t0=0, indf=1)
constraint!(ocp, :initial, [1,0,1], :initial_constraint)
ocots commented 1 week ago
constraint!(ocp, :initial; lb=[1,0,1], ub=[1,0,1], label=:initial_constraint)

Documentation

PierreMartinon commented 1 week ago

Ok, so basically everything after the constraint type is now a named argument. Regarding the doc I would put at least one equality constraint example, currently they all have differnt bounds.

PierreMartinon commented 1 week ago

Another question @ocots :D On the goddard problem, the following box constraint is fine

constraint!(ocp, :state, lb=[r0,v0,mf], ub=[r0+0.2,vmax,m0])

while the functional version fails (restoration error at first iteration)

constraint!(ocp1, :state, f=(x,v)->x, lb=[r0,v0,mf], ub=[r0+0.2,vmax,m0])

I tried splitting in 3 separate constraints, same problem. Did I miss something ?

To clarify: I don't need this specific formulation, I'm just trying to update my test cases with the different constraints types.

jbcaillau commented 1 week ago

@PierreMartinon there is a difference (that should not suffice to alter solving, though) as

Note that while the first syntax is allowed, it is more readable to write

constraint!(ocp, :state; rg=1:end, lb=[r0,v0,mf], ub=[r0+0.2,vmax,m0])

Currently fixing a few things in the constraint! function 🤞🏾

PierreMartinon commented 1 week ago

Yup. As I said, not a blocking issue, I'm just trying to figure out the different possibilities. I'll try to look more closely at the actual constraints in the defined ocp

jbcaillau commented 1 week ago

Yup. As I said, not a blocking issue, I'm just trying to figure out the different possibilities. I'll try to look more closely at the actual constraints in the defined ocp

@PierreMartinon please drop a pointer towards the example not converging with

constraint!(ocp1, :state, f=(x,v)->x, lb=[r0,v0,mf], ub=[r0+0.2,vmax,m0])
PierreMartinon commented 1 week ago

Here you go. The one below (scroll) converges

# general definitions
Cd = 310
Tmax = 3.5
β = 500
b = 2
r0 = 1
v0 = 0
vmax = 0.1
m0 = 1
mf = 0.6
x0 = [ r0, v0, m0 ]
function F0(x)
    r, v, m = x
    D = Cd * v^2 * exp(-β*(r - 1))
    return [ v, -D/m - 1/r^2, 0 ]
end
function F1(x)
    r, v, m = x
    return [ 0, Tmax/m, -b*Tmax ]
end

# box constraints
ocp = Model(variable=true)
state!(ocp, 3)
control!(ocp, 1)
variable!(ocp, 1)
time!(ocp, t0=0, indf=1)
constraint!(ocp, :initial, lb=x0, ub=x0)
constraint!(ocp, :final, rg=3, lb=mf, ub=Inf)
constraint!(ocp, :state, lb=[r0,v0,mf], ub=[r0+0.2,vmax,m0])
constraint!(ocp, :control, lb=0, ub=1)
constraint!(ocp, :variable, lb=0.01, ub=Inf)
objective!(ocp, :mayer, (x0, xf, v) -> xf[1], :max)
dynamics!(ocp, (x, u, v) -> F0(x) + u*F1(x) )

sol = solve(ocp, grid_size=100, print_level=0, tol=1e-8)

but using the functional constraint fails

constraint!(ocp1, :state, f=(x,v)->x, lb=[r0,v0,mf], ub=[r0+0.2,vmax,m0])
PierreMartinon commented 1 week ago

@jbcaillau, with Joseph we may have found the answer: ipopt will internally project the initial guess inside the variable bounds if needed. Here for instance the default init (0.1) will be projected for both r and m.

When using a functional (nonlinear) constraint, there is no such initial guess adjustment, and for this problem in particular, having r<1 (earth radius) is known to be problematic.

I'll recheck with a more sensible initial guess and see if the functional constraint converges.

jbcaillau commented 1 week ago

@PierreMartinon @joseph-gergaud nice! indeed, not having a box prevents Ipopt from making this educated guess.

PierreMartinon commented 6 days ago

Confirmed, we get the convergence back when initializing the state with r>1. No more remarks for me here, you can close when you want.

ocots commented 4 days ago

@jbcaillau Please review what is done: