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)
2 Likes

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.

2 Likes

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.

1 Like

Can you provide your JuMP code?

In addition to working through the logic in maths, I’d encourage you to think about how you can test and validate your code. If you solve the problem and look at the solution, does it respect the constraints? If you change the data, does it work? If you make the red light phase from 0 to tmax is the problem infeasible? And so on.

Sure, I tried to derive the logic again and I find a mistake in setting M, m constants.

Here is updated JuMP code:

using JuMP, Ipopt, Gurobi, Juniper

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

g = 9.81
m = 1
a0 = 1e-2
a1 = 0
a2 = 1.5e-5

# Flat track
alpha = zeros(N+1)

# Dynamics separated per state variable
function f1(x1, x2, u, alpha)
  return 1 / x2
end

function f2(x1, x2, u, alpha)
  return u / (m * x2) - (a0 + a1*x2 + a2*x2^2) / (m * x2) - g * sin(alpha) / x2
end


ipopt = optimizer_with_attributes(Ipopt.Optimizer)
gurobi = optimizer_with_attributes(Gurobi.Optimizer)

model = Model(
    optimizer_with_attributes(
        Juniper.Optimizer,
        "nl_solver" => ipopt,
        "mip_solver" => gurobi,
    ),
)

T = 1400 # Total trip time
v0 = 1 # Initial and minimum speed allowed
vmax = 13.8
eta1 = 1
eta2 = 0

# Mixed-integer related variables
red_phase_start = 400
red_phase_end = 450
tls_index = 200
M1 = red_phase_start
m1 = red_phase_start - T
M2 = T - red_phase_end
m2 = -red_phase_end
epsilon = 1e-4

@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])

# Stations => speed must be small, ideally zero
@constraint(model, x[2, 100] == v0)
@constraint(model, x[1, 100] == 90)
@constraint(model, x[2, 300] == v0)
@constraint(model, x[1, 300] == 325)
@constraint(model, x[2, 400:500] .<= 5) # Speed limit zone
@constraint(model, x[2, 600] == v0)
@constraint(model, x[1, 600] == 770)
@constraint(model, x[2, 900] == v0)
@constraint(model, x[1, 900] == 1040)

# Traffic lights related logic
@variable(model, z[1:2], Bin)
@variable(model, y, Bin) # Encoding time between the start AND the end of the red pahse
@variable(model, delta, Bin) # Binary variable for selecting the vmax

# z1 <=> red_phase_start <= t <= red_phase_end
@constraint(model, red_phase_start - x[1, tls_index] <= (1-z[1])*M1)
@constraint(model, red_phase_start - x[1, tls_index] >= epsilon + (m1 - epsilon)*z[1])
@constraint(model, x[1, tls_index] - red_phase_end <= (1-z[2])*M2)
@constraint(model, x[1, tls_index] - red_phase_end >= epsilon + (m2 - epsilon)*z[2])
# y = z1 AND z2
@constraint(model, y <= z[1])
@constraint(model, y <= z[2])
@constraint(model, y >= z[1] + z[2] - 1)
# if y == 1 then delta = 1 else delta = 0
M = 1
m = 1
@constraint(model, delta <= M*y)
@constraint(model, -delta <= -m*y)
@constraint(model, delta <= 1-m*(1-y))
@constraint(model, -delta <= -1+M*(1-y))
# Choose one of the limits
@constraint(model, x[2, tls_index] <= (1-delta) * vmax + delta * v0)

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

# Dynamics for the first state variable
@NLconstraint(model, [i=1:N-1],
    x[1, i+1] == x[1, i] + 0.5 * h * (
        f1(x[1,i], x[2,i], u[i], alpha[i]) +
        f1(x[1,i+1], x[2,i+1], u[i+1], alpha[i+1])
    )
)

@NLconstraint(model,
    x[1, N+1] == x[1, N] + 0.5 * h * (
        f1(x[1,N], x[2,N], u[N], alpha[N]) +
        f1(x[1,N+1], x[2,N+1], u[N], alpha[N+1])
    )
)

# Dynamics for the second state variable
@NLconstraint(model, [i=1:N-1],
    x[2, i+1] == x[2, i] + 0.5 * h * (
        f2(x[1,i], x[2,i], u[i], alpha[i]) +
        f2(x[1,i+1], x[2,i+1], u[i+1], alpha[i+1])
    )
)

@NLconstraint(model,
    x[2, N+1] == x[2, N] + 0.5 * h * (
        f2(x[1,N], x[2,N], u[N], alpha[N]) +
        f2(x[1,N+1], x[2,N+1], u[N], alpha[N+1])
    )
)

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

λ = 1e-3  # small weight for regularization
@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]) + λ * sum((u[i+1] - u[i])^2 for i = 1:N-1))

# NLP Initialization
set_start_value.(x[1,:], LinRange(0, T, N+1))
set_start_value.(x[2,:], v0*ones(1:N+1))
# Fixing state for stations
fix.(x[:, 1], [0,v0])
fix.(x[:, end], [T, v0])
fix.(x[:, 100], [90, v0])
fix.(x[:, 300], [325, v0])
fix.(x[:, 600], [770, v0])
fix.(x[:, 900], [1040, v0])
optimize!(model)

I must admit that the problem without the mixed-inter part is solved with Ipopt within few iterations, even for a general track (not a flat one).

1 Like