Dependent variable

[FedeClaudi]

Hello, I'm back with some questions.

Assume I have a model with a bunch of variables

    @variables(model, begin
        # CONTROLS
        F_v_max[1] ≤ F_v ≤ F_v_max[2], Infinite(s)        # longitudinal force
        F_ω_max[1] ≤ F_ω ≤ F_ω_max[2], Infinite(s)        # turning force

        # STATE
        -track.width ≤ n ≤ track.width, Infinite(s), (start=0)      # model->track lateral deviation
        -π/2 ≤ ψ ≤ π/2, Infinite(s), (start=0)      # model->track angular deviation

        0 ≤ v ≤ 100, Infinite(s), (start=1)           # speed
        -1e-1 ≤ ω ≤ 1e-1, Infinite(s), (start=0)        # angular velocity

        # time
        0 ≤ t ≤ 60, Infinite(s), (start=0)

        #  kinematics
        θ, Infinite(s)
        x, Infinite(s)
        y, Infinite(s)
        β == 0                                      # slip angle, constant for now

        # constants
        m == 1                                      # mass
    end)

In this model F_v and F_ω are decision variables, in the sense that I want to find the values of these variables that satisfy constraints and minimize an objective function. All the other variables depend on the two F_x controls, and have a specific set of dynamics, which I set with:

    SF = (1 - n * κ(s))/(v * cos(ψ + β) + eps())  # time -> space domain conversion factor

    @constraints(model, begin
        ∂(n, s) == SF * v * sin(ψ + β)
        ∂(ψ, s) == SF * ω - κ(s)
        ∂(v, s) == SF * F_v / m
        # ∂(β, s) == 0             
        ∂(ω, s) == SF * F_ω / m

        ∂(t, s) == SF
        ∂(θ, s) == SF * ω
        ∂(x, s) == SF * v * cos(θ)
        ∂(y, s) == SF * v * sin(θ)
    end)

My understanding is that now all variables with except of F_v and F_ω should not be able to vary independently as their dynamics are set by the derivatives constraints. However I see something weird in the results.
Let's say that this is F_v for the results:

if I do plot(cumsum(value(model[:F_v])) I get:

however this looks very different from plot(value(model[:v])):

Given ∂(v, s) == SF * F_v / m I would expect the two plots to be similar, if not identical. Yet they're very different. Indeed cumsum(value(model[:F_v]) gives values that seem correct to me, while value(model[:v])'s values don't seem to match what you'd expect given F_v. This makes me think that the model is also solving for v, treating it somewhat like a decision variable, instead of simply computing v based on its derivative constraint. Is that the case? Is there a way to have variables in the model, with contraints, but not have the model treat them as decision variables?

Thank you!

[pulsipher]

Welcome back :).

The short answer is this unexpected behavior likely pertains to a modeling error. I do not have all the details of the model, but I can provide some better clarification on what is going on in InfiniteOpt with respect to variables and constraints.

Before I dive into all that, I'll quickly note that β and m are constants so they should not be declared as variables. You can do this but it adds unnecessary complexity for the solver.

First, let me explain that everything defined in the variable macro is treated as a decision variable (the term variable should not be thought of in the classical algebraic sense). When the model is transformed into a finite approximation using a simultaneous direct transcription approach, this becomes a regular NLP that a traditional NLP solver like Ipopt can deal with. Ipopt has no notion of control vs. state variables. It will simply manipulate all the variables simultaneously in accordance with the constraints. The freedom it has to do this is determined by the degrees of freedom inherent in the model formulation. If you have too many degrees of freedom, then the solver will start to manipulate the state variables in ways you probably don't expect. If you have too few degrees of freedom then your problem will likely be infeasible.

InfiniteOpt.jl is a very flexible interface, but the cost of this flexibility is that the responsibility of model accuracy is placed on the user.

Now for your model, here are some thoughts that might be helpful. (Note that I don't have the complete model so I don't have the complete picture).

• First, make sure the problem is solving to optimality, otherwise the output is meaningless.
• If it is optimal, then you probably have too many degrees of freedom so there is probably a modeling error.
• t(s) appears to be unnecessary, you can remove it and its corresponding ODE from the model and compute it after the fact.
• Make sure you have all your initial conditions specified, this is very important. Note that using start does NOT constitute an initial condition, it is an initial guess for the solver.
• Try fixing your controls to a constant value, if the degrees of freedom are correct then this should be a similar response to just simulating the ODEs.
• Double check your math and code to check for modeling bugs.

In the future, I would like to add some automated tools for debugging models like this, but that'll constitute a novel research project in and of itself.

[FedeClaudi]

Thanks for the detailed reply, I really appreciate it.
That clarifies a lot of things actually. I will try to implement all recommendations and let you know.

When you say

make sure the problem is solving to optimality

Does that include increasing the number of iterations such that that's never the cause for ending?
Generally speaking I often end up with a solver has reached a point of infeasibility error, depending on exactly what constraints I put and things like the number of supports.

[pulsipher]

The exact message will depend on the solver. But with Ipopt you want it to say something like "solved to acceptable level" or "optimal". A message that indicates a point of infeasibility means that it wasn't able to find a solution that satisfies all the constraints (i.e., is feasible). This means that the values reported back are essentially garbage.

If you keep getting points of infeasibility this can be caused by two things:

• You have a modeling error.
• The start (i.e., guess) values aren't good enough to get it started.

For the example problem I shared in #219, it got an infeasible solution without the sophisticated guess functions I gave it. This is a common behavior with nonlinear optimization problems that use interior point solvers like Ipopt. One way to get a feasible guess is to simulate your problem in DifferentialEquations.jl with fixed controls, then using that solution as an initial guess.

This is what is motivating the development of an extension package to have InfiniteOpt interface with ModelingToolkit so we can automatically initialize models like this one.

[FedeClaudi]

Hey,

I just wanted to give you an update since you've been helping so much. After much experimentation, what worked was removing those variables that didn't actually contribute to the problem x, y, theta, solve the problem with out and then use DifferentialEquations to forward integrate their ODEs.

This model now converges very nicely:

and navigates around the track

I was hoping to finish the week off with a full solution and the model going around the track... but I think there's something wrong in my ODE integration codes 😄

I think that as fare as InfiniteOpt is concerned, I'm doing everything right now, so I hopefully won't bother you much in the future. Thank you so much for all the help 🙏🏻

[pulsipher]

Awesome, I am glad InfiniteOpt has been helpful. I look forward to seeing your results!