How to incorporate state-dependent state constraint in JuMP

Hi, I am trying to implement energy-efficient optimal control for vehicles, which includes traffic lights in the setup.

The following problem is described as:

\text{minimize}_{u^+, u^-} \int_{0}^{S} \frac{u^+}{\eta_1} - \eta_2u^- \mathrm{d}s
\text{subject to} \quad t'(s) = \frac{1}{v}
\quad \quad \quad \quad \quad v'(s) = \frac{u^+ - u^-}{mv} + \frac{a(s)}{mv},
where s is traveled distance, t is time, v is speed of the vehicle and a(s), where a = g\sin(\alpha(s)) is the gravity-induced acceleration and m is weight of the vehicle.The problem I am facing is that i would like to express constraints on the vehicle speed as a function of distance s (leading to static speed limits based on traveled distance) but also on time t. The dependence on time is motivated by the traffic lights, which are at known distance s with a known red/green phases, something like
v(s) <= V_{max}(s,t(s)) (e.g. at s = 2000 the red phase is for t \in [200, 240], therefore the maximum speed should be zero - or at least small enough).

The problem is discretized using implicit trapezoid method and formulated in JuMP, however I find incorporating the traffic lights constraint challenging. It will probably involve mixed-integer techniques (something like modeling if-else statements). I also tried to wrap it around complementarity constraints, but I guess that won’t fit there.

My question is, how to grasp and formulate the traffic lights constraints in JuMP. Below is MWE (for \alpha = 0).

using JuMP, Ipopt

N = 1000 # Number of discretization points
S = 10e3 # Total distance - 10 km
h = S/N

g = 9.81
m = 1

# Flat track
alpha = zeros(N+1)

# Dynamics
function f(x, u, k)
    return [
      1/x[2],
      u/(m * x[2]) - g*sin(alpha[k])/x[2]
    ]
end

model = Model(Ipopt.Optimizer)
T = 1400 # Total trip time
v0 = 0.5 # Initial & minimum speed
eta1 = 1
eta2 = 0.6

# States & inputs
@variable(model, x[1:2,1:N+1])
@variable(model, u_plus[1:N] >= 0)
@variable(model, u_minus[1:N] >= 0)
@expression(model, u[i=1:N], u_plus[i] - u_minus[i])

# Initial and terminal constraints on state
@constraint(model, x[:,1] == [0, v0])
@constraint(model, x[:, end] == [T, v0])

# Static speed limits
@constraint(model, [i=1:N+1], x[2,i] >= v0)
@constraint(model, [i=1:N+1], x[2,i] <= 9)

# Traction characteristic constraint on input
@constraint(model, [i=1:N], u_plus[i] <= 3/max(x[2,i], 5))
@constraint(model, [i=1:N], u_minus[i] <= 3/max(x[2,i], 5))

# Implicit trapezoid method
@constraint(model, [i=1:N-1], x[:, i+1] == x[:, i] + 1/2 * h * (f(x[:, i], u[i], i) + f(x[:, i+1], u[i+1], i+1)))
@constraint(model, x[:, N+1] == x[:, N] + 1/2 * h *(f(x[:, N], u[N], N) + f(x[:, N+1], u[N], N+1)))

@objective(model, Min, sum([1/2 * h  * ((u_plus[i]/eta1) - (eta2 * u_minus[i]) + (u_plus[i+1]/eta1) - (eta2 * u_minus[i+1])) for i = 1:N-1]))

# Initial guess for NLP
set_start_value.(x[1,:], LinRange(0, T, N+1))
set_start_value.(x[2,:], v0*ones(1:N+1))
fix.(x[:, 1], [0,v0])
fix.(x[:, end], [T, v0])

optimize!(model)
1 Like

I don’t quite understand the “minimize” problem formulas.
What are essentially the decisions? Are they of infinite dimensional (i.e. they are functions)? Or are they ordinary vectors in euclidean n-space? But in the latter case the integrand appears to be pointless. And what does the \prime (t′, v′) mean?

I’m not sure whether it is fitting to cast your optimal control problem as a deterministic optimization problem such that JuMP with solvers can handle.

If you intend to use NLP techniques and aim to attain a suboptimal local solution, then you can use Ipopt. But Ipopt cannot deal with integer decisions.

Moreover, according to my experience, it is not necessarily safe to write logical constraints in a black-box manner. They typically entails big-M, which I assume that solvers will choose it for you (but the quality may depends…).

The decisions are two inputs to the system u^+ and u^-. The meaning of t'(s) and v'(s) is derivative w.r.t. distance, and together, they form a system of differential equations / state-space model (informally called dynamics), which must be satisfied. This continuous time optimal control problem is then discretized into an NLP (the process is called direct transcription).

Hi @stepanoslejsek,

Adding these constraints is going to make it very tricky to solve.

One non-smooth nonlinear approach would be:

s = 2_000
i = round(Int, s / h)
t_s, v_s = x[1,i], x[2,i]
v_max = 100.0
@constraint(model, v_s <= ifelse(200 <= t_s <= 240, 0.0, v_max))

But otherwise, yes, you’ll need to formulate your problem as a MINLP.

cc @pulsipher may have some ideas for how to solve this with InfiniteOpt.jl.

1 Like

Thanks, that is exactly, what I thought. Other way to solve this that I can think of is solving the problem without the constraint and then check, if the time at given position is within the red phase and then enforcing the time to be within previous and/or next green phase. But this branching is going to be needed for every traffic light.

I will try to formulate the problem as MINLP, do you think that Alpine.jl would be sufficient choice?

MINLPs are notoriously difficult to solve. You’ll probably need to try a few different solvers. See also Juniper.jl, EAGO.jl, and if you can get a license, Gurobi.jl.

I tried Juniper.jl, Gurobi.jl and Alpine.jl, all of them failing. Juniper throws an error in the dynamics, where indexing of gradient profile happens → conversion to Int does not solve the problem due to autodiff error. Alpine.jl seems to not support user-defined functions as the error states, that function f cannot be called on variable.

Edit: I managed Juniper.jl to run

I would like to double check, if the MINLP formulation actually makes sense. Let s_i be the position of i-th traffic light, with known red phase from 260 \le t(s_i) \le 300, where t is time.

The bound on the speed v at position s_i is then v(s_i) \le (1 - \delta)v^{\mathrm{max}} + \delta v^{\mathrm{min}}, where \delta \in \{0,1\}.

Now the \delta is enforced to be 1 if t(s_i) is within the red phase interval, otherwise \delta = 0 (IF 260 \le t(s_i) \le 300 THEN \delta = 1 ELSE \delta = 0).

Lets encode condition 260 \le t(s_i) by z_1 and t(s_i) \le 300 by z_2, such that z_1, z_2 \in \{0,1\}. This means that z_1 \iff 260 \le t(s_i) and z_2 \iff t(s_i) \le 300, which can be formulated as:

260 - t(s_i) \le (1-z_1)M_1,
260 - t(s_i) \ge \varepsilon + (m_1 - \varepsilon)z_1,
t(s_i) - 300 \le (1-z_2)M_2,
t(s_i) - 300 \ge \varepsilon + (m_2 - \varepsilon)z_2,

with m_1 = m_2 = -40, M_1 = 260, M_2 = 300 and \varepsilon \approx 10^{-3}. Let now y = z_1 \wedge z_2 formulated as:

y \le z_1,
y \le z_2,
y \ge z_1 + z_2 - 1.

Finally the IF-THEN-ELSE condition is formulated as:

\delta \le My,
-\delta \le -my,
\delta \le 1 - m(1-y),
-\delta \le -1 + M(1-y),

with M = 2 and m = -1.

Is this mixed-integer part formulated correctly? I sincerely appreciate any additional advice, thank you!.

Edit: I realized that z_1 is true only if 260 \le t(s_i), which could be written as t(s_i) \le 260 + Mz_1, same goes for z_2 as -t(s_i) \le -300 + Mz_2 with M = 40.