diff --git a/available_problems_cache.txt b/available_problems_cache.txt new file mode 100644 index 00000000..096f4151 --- /dev/null +++ b/available_problems_cache.txt @@ -0,0 +1,20 @@ +beam +bioreactor +cart_pendulum +chain +dielectrophoretic_particle +double_oscillator +ducted_fan +electric_vehicle +glider +insurance +jackson +moonlander +quadrotor +robbins +robot +rocket +space_shuttle +steering +truck_trailer +vanderpol diff --git a/docs/src/list_of_problems.md b/docs/src/list_of_problems.md index 7cfab0a5..532195d7 100644 --- a/docs/src/list_of_problems.md +++ b/docs/src/list_of_problems.md @@ -33,7 +33,7 @@ The table below summarizes the names and status of the each problem: | dielectrophoretic_particle | ✅ | ✅ | | double_oscillator | ✅ | ✅ | | ducted_fan | ✅ | ✅ | -| electrical_vehicle | ✅ | ✅ | +| electric_vehicle | ✅ | ✅ | | glider | ✅ | ❌ | | insurance | ✅ | ✅ | | jackson | ✅ | ✅ | diff --git a/ext/JuMPModels/beam.jl b/ext/JuMPModels/beam.jl index 9abfb00f..3024fba1 100644 --- a/ext/JuMPModels/beam.jl +++ b/ext/JuMPModels/beam.jl @@ -39,7 +39,7 @@ function OptimalControlProblems.beam(::JuMPBackend; nh::Int=100) end ) - @objective(model, Min, sum(u[t]^2 for t in 0:nh)) + @objective(model, Min, (sum(u[t]^2 for t in 0:nh-1)) / nh) return model end diff --git a/ext/JuMPModels/bioreactor.jl b/ext/JuMPModels/bioreactor.jl index 149b8363..dec18c14 100644 --- a/ext/JuMPModels/bioreactor.jl +++ b/ext/JuMPModels/bioreactor.jl @@ -2,7 +2,7 @@ The Bioreactor Problem: The problem is formulated as a JuMP model and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop) """ -function OptimalControlProblems.bioreactor(::JuMPBackend; nh::Int=100, N::Int=30) +function OptimalControlProblems.bioreactor(::JuMPBackend; nh::Int=250, N::Int=30) # Parameters beta = 1 c = 2 @@ -63,7 +63,7 @@ function OptimalControlProblems.bioreactor(::JuMPBackend; nh::Int=100, N::Int=30 end ) - @objective(model, Max, sum(b[t] / (beta + c) for t in 0:nh)) + @objective(model, Max, step * sum(mu2[t] * b[t] / (beta + c) for t in 0:nh-1)) return model end diff --git a/ext/JuMPModels/cart_pendulum.jl b/ext/JuMPModels/cart_pendulum.jl index f92340c2..bf4930d6 100644 --- a/ext/JuMPModels/cart_pendulum.jl +++ b/ext/JuMPModels/cart_pendulum.jl @@ -4,7 +4,7 @@ The Cart-Pendulum Problem: The objective is to swing the pendulum from the downward position to the upright position in the shortest time possible. The problem is formulated as a JuMP model, and can be found [here](https://arxiv.org/pdf/2303.16746). """ -function OptimalControlProblems.cart_pendulum(::JuMPBackend; nh::Int64=100) +function OptimalControlProblems.cart_pendulum(::JuMPBackend; nh::Int64=250) ## parameters g = 9.81 # gravitation [m/s^2] L = 1.0 # pendulum length [m] @@ -21,13 +21,13 @@ function OptimalControlProblems.cart_pendulum(::JuMPBackend; nh::Int64=100) @variables( model, begin - 0.0 <= tf - ddx - -max_x <= x[0:nh] <= max_x - -max_v <= dx[0:nh] <= max_v - theta[0:nh] - omega[0:nh] - -max_f <= Fex[0:nh] <= max_f + 0.01 <= tf, (start = 0.1) + ddx, (start = 0.1) + -max_x <= x[0:nh] <= max_x, (start = 0.1) + -max_v <= dx[0:nh] <= max_v, (start = 0.1) + theta[0:nh], (start = 0.1) + omega[0:nh], (start = 0.1) + -max_f <= Fex[0:nh] <= max_f, (start = 0.1) end ) @@ -46,14 +46,16 @@ function OptimalControlProblems.cart_pendulum(::JuMPBackend; nh::Int64=100) model, begin step, tf / nh - COG[i=0:nh], L / 2 * [sin(theta[i]) + x[i], -cos(theta[i])] + [x[i], 0] + COG_1[i=0:nh], L / 2 * sin(theta[i]) + x[i] + x[i] + COG_2[i=0:nh], - L / 2 * -cos(theta[i]) alpha[i=0:nh], 1.0 / (I + 0.25 * m * L^2) * 0.5 * L * m * (-ddx * cos(theta[i]) - g * sin(theta[i])) - ddCOG[i=0:nh], L * [-sin(theta[i]), cos(theta[i])] * omega[i] + L / 2 * [cos(theta[i]), sin(theta[i])] * alpha[i] + [ddx, 0] - FXFY[i=0:nh], m * ddCOG[i] + [0, m * g] - eq[i=0:nh], -FXFY[i][1] + Fex[i] - m_cart * ddx - J[i=0:nh], m_cart - c[i=0:nh], eq[i] - J[i] * ddx - ddx_[i=0:nh], -1.0 / J[i] * c[i] + ddCOG_1[i=0:nh], - L * sin(theta[i]) * omega[i] + L / 2 * cos(theta[i]) * alpha[i] + ddx + ddCOG_2[i=0:nh], L * cos(theta[i]) * omega[i] + L / 2 * sin(theta[i]) * alpha[i] + FXFY_1[i=0:nh], m * ddCOG_1[i] + FXFY_2[i=0:nh], m * ddCOG_2[i] + m * g + eq[i=0:nh], -FXFY_1[i] + Fex[i] - m_cart * ddx + c[i=0:nh], eq[i] - m_cart * ddx + ddx_[i=0:nh], -1.0 / m_cart * c[i] alpha_[i=0:nh], 1.0 / (I + 0.25 * m * L^2) * 0.5 * L * m * (-ddx_[i] * cos(theta[i]) - g * sin(theta[i])) end ) diff --git a/ext/JuMPModels/double_oscillator.jl b/ext/JuMPModels/double_oscillator.jl index 8e46bed9..f7661f2d 100644 --- a/ext/JuMPModels/double_oscillator.jl +++ b/ext/JuMPModels/double_oscillator.jl @@ -19,16 +19,16 @@ function OptimalControlProblems.double_oscillator(::JuMPBackend; nh::Int=100) @variables( model, begin - x1[0:nh] - x2[0:nh] - x3[0:nh] - x4[0:nh] - -1.0 <= u[0:nh] <= 1.0 + x1[0:nh], (start = 0.1) + x2[0:nh], (start = 0.1) + x3[0:nh], (start = 0.1) + x4[0:nh], (start = 0.1) + -1.0 <= u[0:nh] <= 1.0, (start = 0.1) end ) # Objective - @objective(model, Min, 0.5 * sum(x1[t]^2 + x2[t]^2 + u[t]^2 for t in 0:nh)) + @objective(model, Min, 0.5 * step * sum(x1[t]^2 + x2[t]^2 + u[t]^2 for t in 0:nh-1)) # Dynamics @expressions( diff --git a/ext/JuMPModels/ducted_fan.jl b/ext/JuMPModels/ducted_fan.jl index 47912754..de36359f 100644 --- a/ext/JuMPModels/ducted_fan.jl +++ b/ext/JuMPModels/ducted_fan.jl @@ -5,7 +5,7 @@ The Ducted Fan Problem: The problem is formulated as a JuMP model. Ref: Graichen, K., & Petit, N. (2009). Incorporating a class of constraints into the dynamics of optimal control problems. Optimal Control Applications and Methods, 30(6), 537-561. """ -function OptimalControlProblems.ducted_fan(::JuMPBackend; nh::Int=100) +function OptimalControlProblems.ducted_fan(::JuMPBackend; nh::Int=250) r = 0.2 # [m] J = 0.05 # [kg.m2] m = 2.2 # [kg] @@ -14,18 +14,16 @@ function OptimalControlProblems.ducted_fan(::JuMPBackend; nh::Int=100) model = Model() - @variable(model, x1[0:nh]) - @variable(model, v1[0:nh]) - @variable(model, x2[0:nh]) - @variable(model, v2[0:nh]) - @variable(model, -deg2rad(30.0) <= α[0:nh] <= deg2rad(30.0)) # radian - @variable(model, vα[0:nh]) - @variable(model, -5.0 <= u1[0:nh] <= 5.0) # [nh] - @variable(model, 0.0 <= u2[0:nh] <= 17.0) # [nh] + @variable(model, x1[0:nh], start = 0.1) + @variable(model, v1[0:nh], start = 0.1) + @variable(model, x2[0:nh], start = 0.1) + @variable(model, v2[0:nh], start = 0.1) + @variable(model, -deg2rad(30.0) <= α[0:nh] <= deg2rad(30.0), start = 0.1) # radian + @variable(model, vα[0:nh], start = 0.1) + @variable(model, -5.0 <= u1[0:nh] <= 5.0, start = 0.1) # [nh] + @variable(model, 0.0 <= u2[0:nh] <= 17.0, start = 0.1) # [nh] @variable(model, 0 <= tf, start = 1.0) - @objective(model, Min, tf / nh * sum(2 * u1[t]^2 + u2[t]^2 for t in 0:nh) + μ * tf) - # Dynamics @expressions( model, @@ -70,5 +68,7 @@ function OptimalControlProblems.ducted_fan(::JuMPBackend; nh::Int=100) end ) + @objective(model, Min, 1. / tf * step * sum(2. * u1[t]^2 + u2[t]^2 for t in 0:nh-1) + μ * tf) + return model end diff --git a/ext/JuMPModels/electrical_vehicle.jl b/ext/JuMPModels/electric_vehicle.jl similarity index 74% rename from ext/JuMPModels/electrical_vehicle.jl rename to ext/JuMPModels/electric_vehicle.jl index 147527f1..08a57482 100644 --- a/ext/JuMPModels/electrical_vehicle.jl +++ b/ext/JuMPModels/electric_vehicle.jl @@ -1,10 +1,10 @@ """ -The Electrical Vehicle Problem - Implement optimal control of an electrical vehicle. +The electric Vehicle Problem + Implement optimal control of an electric vehicle. The problem is formulated as a JuMP model. Ref: [PS2011] Nicolas Petit and Antonio Sciarretta. "Optimal drive of electric vehicles using an inversion-based trajectory generation approach." IFAC Proceedings Volumes 44, no. 1 (2011): 14519-14526. """ -function OptimalControlProblems.electrical_vehicle(::JuMPBackend; nh::Int=100) +function OptimalControlProblems.electric_vehicle(::JuMPBackend; nh::Int=250) D = 10.0 T = 1.0 b1 = 1e3 @@ -18,11 +18,11 @@ function OptimalControlProblems.electrical_vehicle(::JuMPBackend; nh::Int=100) model = Model() - @variable(model, x[0:nh]) - @variable(model, v[0:nh]) - @variable(model, u[0:nh]) + @variable(model, x[0:nh], start = 0.1) + @variable(model, v[0:nh], start = 0.1) + @variable(model, u[0:nh], start = 0.1) - @objective(model, Min, sum(b1 * u[t] * v[t] + b2 * u[t]^2 for t in 0:nh)) + @objective(model, Min, step * sum(b1 * u[t] * v[t] + b2 * u[t]^2 for t in 0:nh-1)) # Dynamics @expressions( diff --git a/ext/JuMPModels/glider.jl b/ext/JuMPModels/glider.jl index bce5482d..f7e31d37 100644 --- a/ext/JuMPModels/glider.jl +++ b/ext/JuMPModels/glider.jl @@ -30,7 +30,7 @@ function OptimalControlProblems.glider(::JuMPBackend; nh::Int64=100) @variables( model, begin - 0 <= t_f, (start = 1.0) + 0 <= tf, (start = 1.0) 0.0 <= x[k=0:nh], (start = x_0 + vx_0 * (k / nh)) y[k=0:nh], (start = y_0 + (k / nh) * (y_f - y_0)) 0.0 <= vx[k=0:nh], (start = vx_0) @@ -44,7 +44,7 @@ function OptimalControlProblems.glider(::JuMPBackend; nh::Int64=100) @expressions( model, begin - step, t_f / nh + step, tf / nh r[i=0:nh], (x[i] / r_0 - 2.5)^2 u[i=0:nh], u_c * (1 - r[i]) * exp(-r[i]) w[i=0:nh], vy[i] - u[i] diff --git a/ext/JuMPModels/insurance.jl b/ext/JuMPModels/insurance.jl index 12c0b5f6..7f6b7a3f 100644 --- a/ext/JuMPModels/insurance.jl +++ b/ext/JuMPModels/insurance.jl @@ -2,7 +2,7 @@ The Insurance Problem: The problem is formulated as a JuMP model and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop) """ -function OptimalControlProblems.insurance(::JuMPBackend; nh::Int=100, N::Int=30) +function OptimalControlProblems.insurance(::JuMPBackend; nh::Int=100) # constants gamma = 0.2 lambda = 0.25 @@ -82,7 +82,7 @@ function OptimalControlProblems.insurance(::JuMPBackend; nh::Int=100, N::Int=30) ) # Objective - @objective(model, Max, sum(U[t] * fx[t] for t in 0:nh)) + @objective(model, Max, step * sum(U[t] * fx[t] for t in 0:nh-1)) return model end diff --git a/ext/JuMPModels/moonlander.jl b/ext/JuMPModels/moonlander.jl index e144cf30..1454221f 100644 --- a/ext/JuMPModels/moonlander.jl +++ b/ext/JuMPModels/moonlander.jl @@ -5,7 +5,7 @@ The Moonlander Problem: The problem is formulated as a JuMP model, and can be found [here](https://arxiv.org/pdf/2303.16746) """ function OptimalControlProblems.moonlander( - ::JuMPBackend; target::Array{Float64}=[5.0, 5.0], nh::Int64=100 + ::JuMPBackend; target::Array{Float64}=[5.0, 5.0], nh::Int64=1000 ) ## parameters if size(target) != (2,) @@ -23,14 +23,14 @@ function OptimalControlProblems.moonlander( @variables( model, begin - 0.0 <= step + 0.01 <= tf, (start = 0.1) # state variables - p1[k=0:nh] - p2[k=0:nh] - dp1[k=0:nh] - dp2[k=0:nh] - theta[k=0:nh] - dtheta[k=0:nh] + p1[k=0:nh], (start = 0.1) + p2[k=0:nh], (start = 0.1) + dp1[k=0:nh], (start = 0.1) + dp2[k=0:nh], (start = 0.1) + theta[k=0:nh], (start = 0.1) + dtheta[k=0:nh], (start = 0.1) # control variables 0 <= F1[k=0:nh] <= max_thrust, (start = 5.0) 0 <= F2[k=0:nh] <= max_thrust, (start = 5.0) @@ -55,6 +55,13 @@ function OptimalControlProblems.moonlander( ) #dynamics + @expressions( + model, + begin + step, tf / nh + end + ) + @expressions( model, begin @@ -95,7 +102,7 @@ function OptimalControlProblems.moonlander( end ) - @objective(model, Min, step * nh) + @objective(model, Min, tf) return model end diff --git a/ext/JuMPModels/quadrotor.jl b/ext/JuMPModels/quadrotor.jl index 1dae5cd0..82595665 100644 --- a/ext/JuMPModels/quadrotor.jl +++ b/ext/JuMPModels/quadrotor.jl @@ -22,17 +22,17 @@ function OptimalControlProblems.quadrotor(::JuMPBackend; nh::Int64=60) @variables( model, begin - 0.0 <= tf - p1[0:nh] - p2[0:nh] - p3[0:nh] - v1[0:nh] - v2[0:nh] - v3[0:nh] + 0.0 <= tf, (start = 0.1) + p1[0:nh], (start = 0.1) + p2[0:nh], (start = 0.1) + p3[0:nh], (start = 0.1) + v1[0:nh], (start = 0.1) + v2[0:nh], (start = 0.1) + v3[0:nh], (start = 0.1) atmin <= at[0:nh] <= atmax, (start = 10.0) - -pi / 2 <= ϕ[0:nh] <= pi / 2 - -pi / 2 <= θ[0:nh] <= pi / 2 - ψ[0:nh] + -pi / 2 <= ϕ[0:nh] <= pi / 2, (start = 0.1) + -pi / 2 <= θ[0:nh] <= pi / 2, (start = 0.1) + ψ[0:nh], (start = 0.1) end ) diff --git a/ext/JuMPModels/robbins.jl b/ext/JuMPModels/robbins.jl index 07d9cafd..36e738c6 100644 --- a/ext/JuMPModels/robbins.jl +++ b/ext/JuMPModels/robbins.jl @@ -2,7 +2,7 @@ The Robbins Problem: The problem is formulated as a JuMP model and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop) """ -function OptimalControlProblems.robbins(::JuMPBackend; nh::Int=100, N::Int=30) +function OptimalControlProblems.robbins(::JuMPBackend; nh::Int=250) # constants alpha = 3 beta = 0 @@ -48,7 +48,7 @@ function OptimalControlProblems.robbins(::JuMPBackend; nh::Int=100, N::Int=30) # Objective @objective( - model, Min, sum(alpha * x1[t] + beta * x1[t]^2 + gamma * u[t]^2 for t in 0:nh) + model, Min, step * sum(alpha * x1[t] + beta * x1[t]^2 + gamma * u[t]^2 for t in 0:nh-1) ) return model diff --git a/ext/JuMPModels/space_Shuttle.jl b/ext/JuMPModels/space_shuttle.jl similarity index 97% rename from ext/JuMPModels/space_Shuttle.jl rename to ext/JuMPModels/space_shuttle.jl index a3791b97..e26a3eab 100644 --- a/ext/JuMPModels/space_Shuttle.jl +++ b/ext/JuMPModels/space_shuttle.jl @@ -1,11 +1,12 @@ """ Space Shuttle Reentry Trajectory Problem: We want to find the optimal trajectory of a space shuttle reentry. - The objective is to minimize the angle of attack at the terminal point. + The objective is to maximize the latitude (cross range) at the terminal point. The problem is formulated as a JuMP model, and can be found [here](https://jump.dev/JuMP.jl/stable/tutorials/nonlinear/space_shuttle_reentry_trajectory/) + Note: no heating limit path constraint """ function OptimalControlProblems.space_shuttle( - ::JuMPBackend; integration_rule::String="rectangular", nh::Int64=503 + ::JuMPBackend; integration_rule::String="trapezoidal", nh::Int64=503 ) ## Global variables w = 203000.0 # weight (lb) diff --git a/ext/JuMPModels/steering.jl b/ext/JuMPModels/steering.jl index f99650af..6a509a92 100644 --- a/ext/JuMPModels/steering.jl +++ b/ext/JuMPModels/steering.jl @@ -23,8 +23,11 @@ function OptimalControlProblems.steering(::JuMPBackend; nh::Int64=100) model = JuMP.Model() @variable(model, u_min <= u[i=1:(nh + 1)] <= u_max, start = 0.0) # control - @variable(model, x[i=1:(nh + 1), j=1:4], start = gen_x0(i, j)) # state - @variable(model, tf, start = 1.0) # final time + @variable(model, x1[i=1:(nh + 1)], start = gen_x0(i, 1)) # state x1 + @variable(model, x2[i=1:(nh + 1)], start = gen_x0(i, 2)) # state x2 + @variable(model, x3[i=1:(nh + 1)], start = gen_x0(i, 3)) # state x3 + @variable(model, x4[i=1:(nh + 1)], start = gen_x0(i, 4)) # state x4 + @variable(model, tf, start = 1.0) # final time @expression(model, h, tf / nh) # step size @objective(model, Min, tf) @@ -35,18 +38,21 @@ function OptimalControlProblems.steering(::JuMPBackend; nh::Int64=100) @constraints( model, begin - con_x1[i=1:nh], x[i + 1, 1] == x[i, 1] + 0.5 * h * (x[i, 3] + x[i + 1, 3]) - con_x2[i=1:nh], x[i + 1, 2] == x[i, 2] + 0.5 * h * (x[i, 4] + x[i + 1, 4]) - con_x3[i=1:nh], - x[i + 1, 3] == x[i, 3] + 0.5 * h * (a * cos(u[i]) + a * cos(u[i + 1])) - con_x4[i=1:nh], - x[i + 1, 4] == x[i, 4] + 0.5 * h * (a * sin(u[i]) + a * sin(u[i + 1])) + con_x1[i=1:nh], x1[i + 1] == x1[i] + 0.5 * h * (x3[i] + x3[i + 1]) + con_x2[i=1:nh], x2[i + 1] == x2[i] + 0.5 * h * (x4[i] + x4[i + 1]) + con_x3[i=1:nh], x3[i + 1] == x3[i] + 0.5 * h * (a * cos(u[i]) + a * cos(u[i + 1])) + con_x4[i=1:nh], x4[i + 1] == x4[i] + 0.5 * h * (a * sin(u[i]) + a * sin(u[i + 1])) end ) # Boundary conditions - @constraint(model, [j = 1:4], x[1, j] == xs[j]) - @constraint(model, [j = 2:4], x[nh + 1, j] == xf[j]) + @constraint(model, x1[1] == xs[1]) + @constraint(model, x2[1] == xs[2]) + @constraint(model, x3[1] == xs[3]) + @constraint(model, x4[1] == xs[4]) + @constraint(model, x2[nh + 1] == xf[2]) + @constraint(model, x3[nh + 1] == xf[3]) + @constraint(model, x4[nh + 1] == xf[4]) return model -end +end \ No newline at end of file diff --git a/ext/JuMPModels/truck_trailer.jl b/ext/JuMPModels/truck_trailer.jl index 73f47b33..3de32997 100644 --- a/ext/JuMPModels/truck_trailer.jl +++ b/ext/JuMPModels/truck_trailer.jl @@ -7,7 +7,7 @@ The Truck Trailer Problem: function OptimalControlProblems.truck_trailer( ::JuMPBackend; data::Array{Float64,2}=[0.4 0.1 0.2; 1.1 0.2 0.2; 0.8 0.1 0.2], - nh::Int64=100, + nh::Int64=200, ) # parameters if size(data) != (3, 3) @@ -40,16 +40,16 @@ function OptimalControlProblems.truck_trailer( model, begin # Final time - 0.0 <= tf + 1.0 <= tf <= 1000, (start = 10) # State variables - x2[0:nh] - y2[0:nh] + x2[0:nh], (start = 0.1) + y2[0:nh], (start = 0.1) -pi / 2 <= theta0[0:nh] <= pi / 2, (start = 0.1) - -pi / 2 <= theta1[0:nh] <= pi / 2, (start = 0.0) - theta2[0:nh], (start = 0.0) + -pi / 2 <= theta1[0:nh] <= pi / 2, (start = 0.1) + theta2[0:nh], (start = 0.1) # Control variables - -0.2 * speedf <= v0[0:nh] <= 0.2 * speedf, (start = -0.2) - -pi / 6 <= delta0[0:nh] <= pi / 6 + -0.2 * speedf <= v0[0:nh] <= 0.2 * speedf, (start = 0.1) + -pi / 6 <= delta0[0:nh] <= pi / 6, (start = 0.1) end ) @@ -70,7 +70,7 @@ function OptimalControlProblems.truck_trailer( begin beta01[j=0:nh], theta0[j] - theta1[j] beta12[j=0:nh], theta1[j] - theta2[j] - step, tf / (nh - 1) + step, tf / nh end ) @constraints( @@ -150,7 +150,7 @@ function OptimalControlProblems.truck_trailer( end ) - @objective(model, Min, tf + sum((beta01[j]^2 + beta12[j]^2) for j in 0:nh)) + @objective(model, Min, tf + step * sum((beta01[j]^2 + beta12[j]^2) for j in 0:nh-1)) return model end diff --git a/ext/JuMPModels/vanderpol.jl b/ext/JuMPModels/vanderpol.jl index e2930445..e5e9e4ca 100644 --- a/ext/JuMPModels/vanderpol.jl +++ b/ext/JuMPModels/vanderpol.jl @@ -13,10 +13,9 @@ function OptimalControlProblems.vanderpol(::JuMPBackend; nh::Int=100) @variables( model, begin - x1[0:nh] - x2[0:nh] - u[0:nh] - step[0:nh] == tf / nh + x1[0:nh], (start = 0.1) + x2[0:nh], (start = 0.1) + u[0:nh], (start = 0.1) end ) @@ -32,6 +31,7 @@ function OptimalControlProblems.vanderpol(::JuMPBackend; nh::Int=100) @expressions( model, begin + step, tf / nh dx1[t=0:nh], x2[t] dx2[t=0:nh], epsilon * omega * (1 - x1[t]^2) * x2[t] - omega^2 * x1[t] + u[t] end @@ -41,12 +41,12 @@ function OptimalControlProblems.vanderpol(::JuMPBackend; nh::Int=100) @constraints( model, begin - con_x1[t=1:nh], x1[t] == x1[t - 1] + 0.5 * step[t] * (dx1[t] + dx1[t - 1]) - con_x2[t=1:nh], x2[t] == x2[t - 1] + 0.5 * step[t] * (dx2[t] + dx2[t - 1]) + con_x1[t=1:nh], x1[t] == x1[t - 1] + 0.5 * step * (dx1[t] + dx1[t - 1]) + con_x2[t=1:nh], x2[t] == x2[t - 1] + 0.5 * step * (dx2[t] + dx2[t - 1]) end ) - @objective(model, Min, sum(0.5 * (x1[t]^2 + x2[t]^2 + u[t]^2) for t in 0:nh)) + @objective(model, Min, step * sum(0.5 * (x1[t]^2 + x2[t]^2 + u[t]^2) for t in 0:nh-1)) return model end diff --git a/ext/MetaData/beam.jl b/ext/MetaData/beam.jl index 276b0c61..8025efff 100644 --- a/ext/MetaData/beam.jl +++ b/ext/MetaData/beam.jl @@ -1 +1,11 @@ -beam_meta = Dict(:name => "beam", :nh => 100, :nvar => 404, :ncon => 305, :minimize => true) +beam_meta = Dict( + :name => "beam", + :nh => 100, + :nvar => 404, + :ncon => 305, + :minimize => true, + :state_name => ["x1", "x2"], + :costate_name => ["con_x1", "con_x2"], + :control_name => ["u"], + :time => ("final_time", "tf", 1.0) +) \ No newline at end of file diff --git a/ext/MetaData/bioreactor.jl b/ext/MetaData/bioreactor.jl index 06d013be..eb16cd79 100644 --- a/ext/MetaData/bioreactor.jl +++ b/ext/MetaData/bioreactor.jl @@ -1,3 +1,11 @@ bioreactor_meta = Dict( - :name => "bioreactor", :nh => 100, :nvar => 505, :ncon => 404, :minimize => false -) + :name => "bioreactor", + :nh => 100, + :nvar => 505, + :ncon => 404, + :minimize => false, + :state_name => ["y", "s", "b"], + :costate_name => ["con_y", "con_s", "con_b"], + :control_name => ["u"], + :time => ("final_time", "T", 300.0) # T = 10 * N where N = 30 by default +) \ No newline at end of file diff --git a/ext/MetaData/cart_pendulum.jl b/ext/MetaData/cart_pendulum.jl index 1f674e45..aae7dd4a 100644 --- a/ext/MetaData/cart_pendulum.jl +++ b/ext/MetaData/cart_pendulum.jl @@ -1,3 +1,11 @@ cart_pendulum_meta = Dict( - :name => "cart_pendulum", :nh => 100, :nvar => 507, :ncon => 405, :minimize => true + :name => "cart_pendulum", + :nh => 100, + :nvar => 507, + :ncon => 405, + :minimize => true, + :state_name => ["x", "dx", "theta", "omega"], + :costate_name => ["d_x", "d_dx", "d_theta", "d_omega"], + :control_name => ["Fex"], + :time => ("final_time", "tf", nothing) ) diff --git a/ext/MetaData/chain.jl b/ext/MetaData/chain.jl index f09be65c..8e93ae56 100644 --- a/ext/MetaData/chain.jl +++ b/ext/MetaData/chain.jl @@ -1,3 +1,11 @@ chain_meta = Dict( - :name => "chain", :nh => 100, :nvar => 404, :ncon => 305, :minimize => true -) + :name => "chain", + :nh => 100, + :nvar => 404, + :ncon => 305, + :minimize => true, + :state_name => ["x1", "x2", "x3"], + :costate_name => ["con_x1", "con_x2", "con_x3"], + :control_name => ["u"], + :time => ("final_time", "tf", 1.0) +) \ No newline at end of file diff --git a/ext/MetaData/dielectrophoretic_particle.jl b/ext/MetaData/dielectrophoretic_particle.jl index b0552c38..0084aedd 100644 --- a/ext/MetaData/dielectrophoretic_particle.jl +++ b/ext/MetaData/dielectrophoretic_particle.jl @@ -4,4 +4,8 @@ dielectrophoretic_particle_meta = Dict( :nvar => 304, :ncon => 203, :minimize => true, -) + :state_name => ["x", "y"], + :costate_name => ["con_x", "con_y"], + :control_name => "u", + :time => ("final_time", "tf", nothing) +) \ No newline at end of file diff --git a/ext/MetaData/double_oscillator.jl b/ext/MetaData/double_oscillator.jl index 67b37d9a..d413fd9e 100644 --- a/ext/MetaData/double_oscillator.jl +++ b/ext/MetaData/double_oscillator.jl @@ -1,7 +1,11 @@ double_oscillator_meta = Dict( :name => "double_oscillator", - :nh => nothing, + :nh => 100, #:nh => nothing, :nvar => nothing, :ncon => nothing, :minimize => true, -) + :state_name => ["x1", "x2", "x3", "x4"], + :costate_name => ["con_x1", "con_x2", "con_x3", "con_x4"], + :control_name => "u", + :time => ("final_time", "tf", 2π) +) \ No newline at end of file diff --git a/ext/MetaData/ducted_fan.jl b/ext/MetaData/ducted_fan.jl index b7022db6..de038e1e 100644 --- a/ext/MetaData/ducted_fan.jl +++ b/ext/MetaData/ducted_fan.jl @@ -1,7 +1,11 @@ ducted_fan_meta = Dict( :name => "ducted_fan", - :nh => nothing, + :nh => 100, # :nh => nothing, :nvar => nothing, :ncon => nothing, :minimize => true, -) + :state_name => ["x1", "v1", "x2", "v2", "α", "vα"], + :costate_name => ["con_x1", "con_v1", "con_x2", "con_v2", "con_α", "con_vα"], + :control_name => ["u1", "u2"], + :time => ("final_time", "tf", nothing) +) \ No newline at end of file diff --git a/ext/MetaData/electric_vehicle.jl b/ext/MetaData/electric_vehicle.jl new file mode 100644 index 00000000..ea3a78b3 --- /dev/null +++ b/ext/MetaData/electric_vehicle.jl @@ -0,0 +1,11 @@ +electric_vehicle_meta = Dict( + :name => "electric_vehicle", + :nh => 100, #:nh => nothing, + :nvar => nothing, + :ncon => nothing, + :minimize => true, + :state_name => ["x", "v"], + :costate_name => ["cond_x", "cond_v"], + :control_name => "u", + :time => ("final_time", "tf", 1.0) +) \ No newline at end of file diff --git a/ext/MetaData/electrical_vehicle.jl b/ext/MetaData/electrical_vehicle.jl deleted file mode 100644 index 50085eea..00000000 --- a/ext/MetaData/electrical_vehicle.jl +++ /dev/null @@ -1,7 +0,0 @@ -electrical_vehicle_meta = Dict( - :name => "electrical_vehicle", - :nh => nothing, - :nvar => nothing, - :ncon => nothing, - :minimize => true, -) diff --git a/ext/MetaData/glider.jl b/ext/MetaData/glider.jl index 2495ebbd..f956c8e8 100644 --- a/ext/MetaData/glider.jl +++ b/ext/MetaData/glider.jl @@ -1,3 +1,11 @@ glider_meta = Dict( - :name => "glider", :nh => 100, :nvar => 506, :ncon => 407, :minimize => false + :name => "glider", + :nh => 100, + :nvar => 506, + :ncon => 407, + :minimize => false, + :state_name => ["x", "y", "vx", "vy"], + :costate_name => ["x_eqn", "y_eqn", "vx_eqn", "vy_eqn"], + :control_name => ["cL"], + :time => ("final_time", "tf", nothing) ) diff --git a/ext/MetaData/insurance.jl b/ext/MetaData/insurance.jl index 1f0084ec..87c14852 100644 --- a/ext/MetaData/insurance.jl +++ b/ext/MetaData/insurance.jl @@ -1,3 +1,11 @@ insurance_meta = Dict( - :name => "insurance", :nh => 100, :nvar => 910, :ncon => 809, :minimize => false -) + :name => "insurance", + :nh => 100, #:nh => nothing, + :nvar => nothing, + :ncon => nothing, + :minimize => false, + :state_name => ["I", "m", "x3"], + :costate_name => ["con_dI", "con_dm", "con_dx3"], + :control_name => ["h", "R", "H", "U", "dUdR"], + :time => ("final_time", "tf", 10) +) \ No newline at end of file diff --git a/ext/MetaData/jackson.jl b/ext/MetaData/jackson.jl index bde0add9..32302f6d 100644 --- a/ext/MetaData/jackson.jl +++ b/ext/MetaData/jackson.jl @@ -1,3 +1,11 @@ jackson_meta = Dict( - :name => "jackson", :nh => 100, :nvar => 404, :ncon => 303, :minimize => false -) + :name => "jackson", + :nh => 100, + :nvar => 404, + :ncon => 303, + :minimize => false, + :state_name => ["a", "b", "x3"], + :costate_name => ["con_da", "con_db", "con_dx3"], + :control_name => ["u"], + :time => ("final_time", "tf", 4.0) +) \ No newline at end of file diff --git a/ext/MetaData/moonlander.jl b/ext/MetaData/moonlander.jl index 0a17cdcf..d5266889 100644 --- a/ext/MetaData/moonlander.jl +++ b/ext/MetaData/moonlander.jl @@ -1,3 +1,11 @@ moonlander_meta = Dict( - :name => "moonlander", :nh => 100, :nvar => 610, :ncon => 809, :minimize => true -) + :name => "moonlander", + :nh => 100, + :nvar => 610, + :ncon => 809, + :minimize => true, + :state_name => ["p1", "p2", "dp1", "dp2", "theta", "dtheta"], + :costate_name => ["d_p1", "d_p2", "d_dp1", "d_dp2", "d_theta", "d_dtheta"], + :control_name => ["F1", "F2"], + :time => ("final_time", "tf", nothing) +) \ No newline at end of file diff --git a/ext/MetaData/quadrotor.jl b/ext/MetaData/quadrotor.jl index 2da328b9..5abf6c2d 100644 --- a/ext/MetaData/quadrotor.jl +++ b/ext/MetaData/quadrotor.jl @@ -1,7 +1,11 @@ quadrotor_meta = Dict( :name => "quadrotor", - :nh => nothing, + :nh => 60, #:nh => nothing, :nvar => nothing, :ncon => nothing, :minimize => true, -) + :state_name => ["p1", "p2", "p3", "v1", "v2", "v3"], + :costate_name => ["d_p1", "d_p2", "d_p3", "d_v1", "d_v2", "d_v3"], + :control_name => ["at", "ϕ", "θ", "ψ"], + :time => ("final_time", "tf", nothing) +) \ No newline at end of file diff --git a/ext/MetaData/robbins.jl b/ext/MetaData/robbins.jl index ddb9fbd8..b5332bd6 100644 --- a/ext/MetaData/robbins.jl +++ b/ext/MetaData/robbins.jl @@ -1,3 +1,11 @@ robbins_meta = Dict( - :name => "robbins", :nh => 100, :nvar => 505, :ncon => 407, :minimize => true -) + :name => "robbins", + :nh => 100, + :nvar => 505, + :ncon => 407, + :minimize => true, + :state_name => ["x1", "x2", "x3"], + :costate_name => ["con_x1", "con_x2", "con_dx3"], + :control_name => ["u"], + :time => ("final_time", "tf", 10) +) \ No newline at end of file diff --git a/ext/MetaData/robot.jl b/ext/MetaData/robot.jl index 653962dd..070488f9 100644 --- a/ext/MetaData/robot.jl +++ b/ext/MetaData/robot.jl @@ -1,3 +1,11 @@ robot_meta = Dict( - :name => "robot", :nh => 100, :nvar => 910, :ncon => 612, :minimize => true -) + :name => "robot", + :nh => 100, + :nvar => 910, + :ncon => 612, + :minimize => true, + :state_name => ["rho", "rho_dot", "the", "the_dot", "phi", "phi_dot"], + :costate_name => ["con_rho", "con_rho_dot", "con_the", "con_the_dot", "con_phi", "con_phi_dot"], + :control_name => ["u_rho", "u_the", "u_phi"], + :time => ("final_time", "tf", nothing) +) \ No newline at end of file diff --git a/ext/MetaData/rocket.jl b/ext/MetaData/rocket.jl index 0ec4c5bd..243214c8 100644 --- a/ext/MetaData/rocket.jl +++ b/ext/MetaData/rocket.jl @@ -1,3 +1,11 @@ rocket_meta = Dict( - :name => "rocket", :nh => 100, :nvar => 405, :ncon => 304, :minimize => false -) + :name => "rocket", + :nh => 100, + :nvar => 405, + :ncon => 304, + :minimize => false, + :state_name => ["h", "v", "m"], + :costate_name => ["con_dh", "con_dv", "con_dm"], + :control_name => ["T"], + :time => ("step", "step", nothing) +) \ No newline at end of file diff --git a/ext/MetaData/space_shuttle.jl b/ext/MetaData/space_shuttle.jl index 6bf6d842..5689ca8d 100644 --- a/ext/MetaData/space_shuttle.jl +++ b/ext/MetaData/space_shuttle.jl @@ -1,7 +1,11 @@ space_shuttle_meta = Dict( :name => "space_shuttle", - :nh => nothing, + :nh => 503, #:nh => nothing, :nvar => nothing, :ncon => nothing, :minimize => false, -) + :state_name => ["scaled_h", "ϕ", "θ", "scaled_v", "γ", "ψ"], + :costate_name => ["con_dh", "con_dϕ", "con_dθ", "con_dv", "con_dγ", "con_dψ"], + :control_name => ["α", "β"], + :time => ("step", "Δt", nothing) +) \ No newline at end of file diff --git a/ext/MetaData/steering.jl b/ext/MetaData/steering.jl index afb1fe98..8eeb0491 100644 --- a/ext/MetaData/steering.jl +++ b/ext/MetaData/steering.jl @@ -1,3 +1,11 @@ steering_meta = Dict( - :name => "steering", :nh => 100, :nvar => 506, :ncon => 408, :minimize => true -) + :name => "steering", + :nh => 100, + :nvar => 506, + :ncon => 408, + :minimize => true, + :state_name => ["x1", "x2", "x3", "x4"], + :costate_name => ["con_x1", "con_x2", "con_x3", "con_x4"], + :control_name => "u", + :time => ("final_time", "tf", nothing) +) \ No newline at end of file diff --git a/ext/MetaData/truck_trailer.jl b/ext/MetaData/truck_trailer.jl index 79f4e134..f0f16ece 100644 --- a/ext/MetaData/truck_trailer.jl +++ b/ext/MetaData/truck_trailer.jl @@ -1,3 +1,11 @@ truck_trailer_meta = Dict( - :name => "truck", :nh => nothing, :nvar => nothing, :ncon => nothing, :minimize => true -) + :name => "truck_trailer", + :nh => 100, #:nh => nothing, + :nvar => nothing, + :ncon => nothing, + :minimize => true, + :state_name => ["x2", "y2", "theta2", "theta1", "theta0"], + :costate_name => ["d_x2", "d_y2", "d_theta2", "d_theta1", "d_theta0"], + :control_name => ["v0", "delta0"], + :time => ("final_time", "tf", nothing) +) \ No newline at end of file diff --git a/ext/MetaData/vanderpol.jl b/ext/MetaData/vanderpol.jl index 47d6bb48..d216423e 100644 --- a/ext/MetaData/vanderpol.jl +++ b/ext/MetaData/vanderpol.jl @@ -1,3 +1,11 @@ vanderpol_meta = Dict( - :name => "vanderpol", :nh => 100, :nvar => 404, :ncon => 303, :minimize => true -) + :name => "vanderpol", + :nh => 100, + :nvar => 404, + :ncon => 303, + :minimize => true, + :state_name => ["x1", "x2"], + :costate_name => ["con_x1", "con_x2"], + :control_name => "u", + :time => ("final_time", "tf", 2.0) +) \ No newline at end of file diff --git a/ext/OptimalControlModels/beam.jl b/ext/OptimalControlModels/beam.jl index 47d661e2..22f97c13 100644 --- a/ext/OptimalControlModels/beam.jl +++ b/ext/OptimalControlModels/beam.jl @@ -20,7 +20,5 @@ function OptimalControlProblems.beam(::OptimalControlBackend; nh::Int=100) init = (state=[0.0, 0.0], control=0.0) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/bioreactor.jl b/ext/OptimalControlModels/bioreactor.jl index 4c3bc907..d598f416 100644 --- a/ext/OptimalControlModels/bioreactor.jl +++ b/ext/OptimalControlModels/bioreactor.jl @@ -2,17 +2,26 @@ The Bioreactor Problem: The problem is formulated as an OptimalControl model and can be found [here](https://github.com/control-toolbox/bocop/tree/main/bocop) """ -function OptimalControlProblems.bioreactor(::OptimalControlBackend; nh::Int=100, N::Int=30) - # constants - beta = 1 - c = 2 - gamma = 1 - halfperiod = 5 - Ks = 0.05 - mu2m = 0.1 - mubar = 1 - r = 0.005 - T = 10 * N +function OptimalControlProblems.bioreactor(::OptimalControlBackend; nh::Int=250, N::Int=30) + + # METHANE PROBLEM + # mu2 according to growth model + # mu according to light model + # time scale is [0,10] for 24h (day then night) + + # growth model MONOD + function growth(s, mu2m, Ks) + return mu2m * s / (s + Ks) + end + + # light model: max^2 (0,sin) * mubar + # DAY/NIGHT CYCLE: [0,2 halfperiod] rescaled to [0,2pi] + function light(time, halfperiod) + pi = 3.141592653589793 + days = time / (halfperiod * 2) + tau = (days - floor(days)) * 2 * pi + return max(0, sin(tau))^2 + end # Model ocp = @def begin @@ -34,41 +43,24 @@ function OptimalControlProblems.bioreactor(::OptimalControlBackend; nh::Int=100, y = x[1] s = x[2] b = x[3] - mu2 = mu2m * s(t) / (s(t) + Ks) + mu = light(t, halfperiod) * mubar + mu2 = growth(s(t), mu2m, Ks) [0, 0, 0.001] ≤ x(t) ≤ [Inf, Inf, Inf] 0 ≤ u(t) ≤ 1 0.05 ≤ y(0) ≤ 0.25 0.5 ≤ s(0) ≤ 5 0.5 ≤ b(0) ≤ 3 - - # dynamics - ẋ(t) == dynamics(t, x(t), u(t)) - - # objective - ∫(b(t) / (beta + c)) → max - end - - # dynamics - function dynamics(t, x, u) - y, s, b = x - pi = 3.141592653589793 - days = t / (halfperiod * 2) - tau = (days - floor(days)) * 2 * pi - light = max(0, sin(tau))^2 - mu = light * mubar - mu2 = mu2m * s / (s + Ks) - return [ - mu * y / (1 + y) - (r + u) * y, - -mu2 * b + u * beta * (gamma * y - s), - (mu2 - u * beta) * b, + ẋ(t) == [ + mu * y(t) / (1 + y(t)) - (r + u(t)) * y(t), + -mu2 * b(t) + u(t) * beta * (gamma * y(t) - s(t)), + (mu2 - u(t) * beta) * b(t), ] + ∫(mu2 * b(t) / (beta + c)) → max end # Initial guess init = (state=[50, 50, 50], control=0.5) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/cart_pendulum.jl b/ext/OptimalControlModels/cart_pendulum.jl index 8c0ab556..50ded98a 100644 --- a/ext/OptimalControlModels/cart_pendulum.jl +++ b/ext/OptimalControlModels/cart_pendulum.jl @@ -90,10 +90,11 @@ function OptimalControlProblems.cart_pendulum(::OptimalControlBackend; nh::Int=1 end # initial guess - init = () + xinit = t -> [0.1, 0.1, 0.1, 0.1] # [x1, dx, theta, omega] + uinit = [0.1] # [Fex] + varinit = [0.1, 0.1] # [tf, ddx] + init = (state=xinit, control=uinit, variable=varinit) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/chain.jl b/ext/OptimalControlModels/chain.jl index 440de72d..e92defa0 100644 --- a/ext/OptimalControlModels/chain.jl +++ b/ext/OptimalControlModels/chain.jl @@ -58,7 +58,5 @@ function OptimalControlProblems.chain(::OptimalControlBackend; nh::Int=100) init = (state=xinit, control=uinit) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/dielectrophoretic_particle.jl b/ext/OptimalControlModels/dielectrophoretic_particle.jl index 9a0a1309..70a08335 100644 --- a/ext/OptimalControlModels/dielectrophoretic_particle.jl +++ b/ext/OptimalControlModels/dielectrophoretic_particle.jl @@ -62,7 +62,5 @@ function OptimalControlProblems.dielectrophoretic_particle( init = (state=[1.0, 1.0], control=0.0, variable=1.0) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/double_oscillator.jl b/ext/OptimalControlModels/double_oscillator.jl index 8db56f32..b6fdacb9 100644 --- a/ext/OptimalControlModels/double_oscillator.jl +++ b/ext/OptimalControlModels/double_oscillator.jl @@ -55,10 +55,10 @@ function OptimalControlProblems.double_oscillator(::OptimalControlBackend; nh::I end # Initial guess - init = () + xinit = t -> [0.1, 0.1, 0.1, 0.1] # [x1, x2, x3, x4] + uinit = [0.1] # [u] + init = (state=xinit, control=uinit) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/ducted_fan.jl b/ext/OptimalControlModels/ducted_fan.jl index 067022e9..6b393b81 100644 --- a/ext/OptimalControlModels/ducted_fan.jl +++ b/ext/OptimalControlModels/ducted_fan.jl @@ -77,10 +77,11 @@ function OptimalControlProblems.ducted_fan(::OptimalControlBackend; nh::Int=100) end # Initial guess - init = () + xinit = t -> [0.1, 0.1, 0.1, 0.1, 0.1, 0.1] # [x1, v1, x2, v2, α, vα] + uinit = [0.1, 0.1] # [u1, u2] + varinit = [1.0] # [tf] + init = (state=xinit, control=uinit, variable=varinit) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/electrical_vehicle.jl b/ext/OptimalControlModels/electric_vehicle.jl similarity index 74% rename from ext/OptimalControlModels/electrical_vehicle.jl rename to ext/OptimalControlModels/electric_vehicle.jl index 0941e203..d34cadb2 100644 --- a/ext/OptimalControlModels/electrical_vehicle.jl +++ b/ext/OptimalControlModels/electric_vehicle.jl @@ -1,10 +1,10 @@ """ -The Electrical Vehicle Problem - Implement optimal control of an electrical vehicle. +The electric Vehicle Problem + Implement optimal control of an electric vehicle. The problem is formulated as an OptimalControl model. Ref: [PS2011] Nicolas Petit and Antonio Sciarretta. "Optimal drive of electric vehicles using an inversion-based trajectory generation approach." IFAC Proceedings Volumes 44, no. 1 (2011): 14519-14526. """ -function OptimalControlProblems.electrical_vehicle(::OptimalControlBackend; nh::Int=100) +function OptimalControlProblems.electric_vehicle(::OptimalControlBackend; nh=250) # parameters D = 10.0 tf = 1.0 @@ -16,14 +16,7 @@ function OptimalControlProblems.electrical_vehicle(::OptimalControlBackend; nh:: p0, p1, p2, p3 = (3.0, 0.4, -1.0, 0.1) ocp = @def begin - ## parameters - D = 10.0 - tf = 1.0 - b1 = 1e3 - b2 = 1e3 - h0 = 0.1 - h1 = 1.0 - h2 = 1e-3 + ## define the problem t ∈ [0.0, tf], time x ∈ R², state @@ -59,10 +52,10 @@ function OptimalControlProblems.electrical_vehicle(::OptimalControlBackend; nh:: end ## Initial guess - init = () + xinit = t -> [0.1, 0.1] # [pos, v] + uinit = [0.1] # [u] + init = (state=xinit, control=uinit) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/glider.jl b/ext/OptimalControlModels/glider.jl index 194f5105..ebbeafe3 100644 --- a/ext/OptimalControlModels/glider.jl +++ b/ext/OptimalControlModels/glider.jl @@ -3,6 +3,7 @@ Hang Glider Problem: We want to find the optimal trajectory of a hang glider. The objective is to maximize the final horizontal position of the glider while in the presence of a thermal updraft. The problem is formulated as an OptimalControl model. + Original formulation from MadNLP/COPSBenchmark """ function OptimalControlProblems.glider(::OptimalControlBackend; nh::Int=100) # parameters @@ -25,31 +26,28 @@ function OptimalControlProblems.glider(::OptimalControlBackend; nh::Int=100) cL_max = 1.4 t0 = 0.0 + function dynamics(x, u) + x1, y, vx, vy = x + cL = u + ## Helper functions + r = (x1 / r_0 - 2.5)^2 + UpD = u_c * (1 - r) * exp(-r) + w = vy - UpD + v = sqrt(vx^2 + w^2) + D = 0.5 * (c0 + c1 * (cL^2)) * rho * S * (v^2) + L = 0.5 * cL * rho * S * (v^2) + vx_dot = (-L * (w / v) - D * (vx / v)) / m + vy_dot = ((L * (vx / v) - D * (w / v)) / m) - g + return [vx, vy, vx_dot, vy_dot] + end + ocp = @def begin - ## parameters - x_0 = 0.0 - y_0 = 1000.0 - y_f = 900.0 - vx_0 = 13.23 - vx_f = 13.23 - vy_0 = -1.288 - vy_f = -1.288 - u_c = 2.5 - r_0 = 100.0 - m = 100.0 - g = 9.81 - c0 = 0.034 - c1 = 0.069662 - S = 14.0 - rho = 1.13 - cL_min = 0.0 - cL_max = 1.4 - t0 = 0.0 + ## define the problem - tf ∈ R¹, variable + tf ∈ R, variable t ∈ [t0, tf], time x ∈ R⁴, state - u ∈ R¹, control + u ∈ R, control ## state variables y = x₂ vx = x₃ @@ -69,6 +67,7 @@ function OptimalControlProblems.glider(::OptimalControlBackend; nh::Int=100) vx(t0) == vx_0, (vx0_con) vy(t0) == vy_0, (vy0_con) # final conditions + tf ≥ 10.0 y(tf) == y_f, (yf_con) vx(tf) == vx_f, (vxf_con) vy(tf) == vy_f, (vyf_con) @@ -80,32 +79,12 @@ function OptimalControlProblems.glider(::OptimalControlBackend; nh::Int=100) x₁(tf) → max end - function dynamics(x, u) - x1, y, vx, vy = x - cL = u - ## Helper functions - r = (x1 / r_0 - 2.5)^2 - UpD = u_c * (1 - r) * exp(-r) - w = vy - UpD - v = sqrt(vx^2 + w^2) - D = 0.5 * (c0 + c1 * (cL^2)) * rho * S * (v^2) - L = 0.5 * cL * rho * S * (v^2) - vx_dot = (-L * (w / v) - D * (vx / v)) / m - vy_dot = ((L * (vx / v) - D * (w / v)) / m) - g - return [vx, vy, vx_dot, vy_dot] - end - # Initial guess - tf = 90.0 - xinit = [ - [x_0 + vx_0 * (k / nh), y_0 + (k / nh) * (y_f - y_0), vx_0, vy_0] for k in 0:nh - ] + tf = (y_f - y_0) / vy_0 + xinit = t -> [x_0 + vx_0 * t, y_0 + t / tf * (y_f - y_0), vx_0, vy_0] uinit = cL_max / 2.0 - time_vec = LinRange(0.0, tf, nh + 1) - init = (time=time_vec, state=xinit, control=uinit, variable=1.0) + init = (state=xinit, control=uinit, variable=tf) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/insurance.jl b/ext/OptimalControlModels/insurance.jl index 62967937..0c939eeb 100644 --- a/ext/OptimalControlModels/insurance.jl +++ b/ext/OptimalControlModels/insurance.jl @@ -71,10 +71,11 @@ function OptimalControlProblems.insurance(::OptimalControlBackend; nh::Int=100) end # Initial guess - init = () + xinit = t -> [0.1, 0.1, 0.1] # [I, m, x3] + uinit = [0.1, 0.1, 0.1, 0.1, 0.1] # [h, R, H, U, dUdR] + varinit = [0.1] # [P] + init = (state=xinit, control=uinit, variable=varinit) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/jackson.jl b/ext/OptimalControlModels/jackson.jl index b11bafa9..02e332cc 100644 --- a/ext/OptimalControlModels/jackson.jl +++ b/ext/OptimalControlModels/jackson.jl @@ -27,10 +27,10 @@ function OptimalControlProblems.jackson(::OptimalControlBackend; nh::Int=100, N: end # Initial guess - init = () + xinit = t -> [0.1, 0.1, 0.1] # [a, b, x3] + uinit = [0.1] # [u] + init = (state=xinit, control=uinit) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/moonlander.jl b/ext/OptimalControlModels/moonlander.jl index f3b90efc..765af53d 100644 --- a/ext/OptimalControlModels/moonlander.jl +++ b/ext/OptimalControlModels/moonlander.jl @@ -5,25 +5,34 @@ The Moonlander Problem: The problem is formulated as an OptimalControl model. """ function OptimalControlProblems.moonlander( - ::OptimalControlBackend; target::Array{Float64}=[5.0, 5.0], nh::Int=100 + ::OptimalControlBackend; p_f::Array{Float64}=[5.0, 5.0], nh::Int=100 ) ## parameters - if size(target) != (2,) - error("The input matrix must be 3x3.") - end m = 1.0 g = 9.81 I = 0.1 D = 1.0 max_thrust = 2 * g + # dynamics + function dynamics(x, u) + p1, p2, dp1, dp2, theta, dtheta = x + F1, F2 = u + + F_r = [ + cos(theta) -sin(theta) p1 + sin(theta) cos(theta) p2 + 0.0 0.0 1.0 + ] + F_tot = (F_r * [0; F1 + F2; 0])[1:2] + ddp1 = (1 / m) * F_tot[1] + ddp2 = (1 / m) * F_tot[2] - g + ddtheta = (1 / I) * (D / 2) * (F2 - F1) + + return [dp1, dp2, ddp1, ddp2, dtheta, ddtheta] + end + ocp = @def begin - ## parameters - m = 1.0 - g = 9.81 - I = 0.1 - D = 1.0 - max_thrust = 2 * g ## define the problem tf ∈ R, variable @@ -55,8 +64,9 @@ function OptimalControlProblems.moonlander( theta(0.0) == 0.0, (theta_ic) dtheta(0.0) == 0.0, (dtheta_ic) # final conditions - p1(tf) == target[1], (p1_fc) - p2(tf) == target[2], (p2_fc) + tf >= 0.1 + p1(tf) == p_f[1], (p1_fc) + p2(tf) == p_f[2], (p2_fc) dp1(tf) == 0.0, (dp1_fc) dp2(tf) == 0.0, (dp2_fc) @@ -67,29 +77,12 @@ function OptimalControlProblems.moonlander( tf → min end - # dynamics - function dynamics(x, u) - p1, p2, dp1, dp2, theta, dtheta = x - F1, F2 = u - - F_r = [ - cos(theta) -sin(theta) p1 - sin(theta) cos(theta) p2 - 0.0 0.0 1.0 - ] - F_tot = (F_r * [0; F1 + F2; 0])[1:2] - ddp1 = (1 / m) * F_tot[1] - ddp2 = (1 / m) * F_tot[2] - g - ddtheta = (1 / I) * (D / 2) * (F2 - F1) - - return [dp1, dp2, ddp1, ddp2, dtheta, ddtheta] - end - # Initial guess - init = (control=[5.0, 5.0],) + xinit = t -> [0.1, 0.1, 0.1, 0.1, 0.1, 0.1] # [p1, p2, dp1, dp2, theta, dtheta] + uinit = [5.0, 5.0] # [F1, F2] + varinit = [0.1] # [tf] + init = (state=xinit, control=uinit, variable=varinit) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/quadrotor.jl b/ext/OptimalControlModels/quadrotor.jl index 111e4662..461b3c75 100644 --- a/ext/OptimalControlModels/quadrotor.jl +++ b/ext/OptimalControlModels/quadrotor.jl @@ -17,21 +17,32 @@ function OptimalControlProblems.quadrotor(::OptimalControlBackend; nh::Int=60) pf = [0.01, 5.0, 2.5] vf = [0.0, 0.0, 0.0] + function dynamics(x, u) + p1, p2, p3, v1, v2, v3, ϕ, θ = x + at, ϕ_dot, θ_dot, ψ = u + + cr = cos(ϕ) + sr = sin(ϕ) + cp = cos(θ) + sp = sin(θ) + cy = cos(ψ) + sy = sin(ψ) + R = [ + (cy*cp) (cy * sp * sr-sy * cr) (cy * sp * cr+sy * sr) + (sy*cp) (sy * sp * sr+cy * cr) (sy * sp * cr-cy * sr) + (-sp) (cp*sr) (cp*cr) + ] + at_ = R * [0; 0; at] + g_ = [0; 0; -g] + a = g_ + at_ + + return [v1, v2, v3, a[1], a[2], a[3], ϕ_dot, θ_dot] + end + ocp = @def begin - ## parameters - g = 9.81 - atmin = 0 - atmax = 9.18 * 5 - tiltmax = 1.1 / 2 - dtiltmax = 6.0 / 2 - p0 = [0.0, 0.0, 2.5] - v0 = [0, 0, 0] - u0 = [9.81, 0, 0, 0] - pf = [0.01, 5.0, 2.5] - vf = [0.0, 0.0, 0.0] ## define the problem - tf ∈ R¹, variable + tf ∈ R, variable t ∈ [0.0, tf], time x ∈ R⁸, state u ∈ R⁴, control @@ -53,14 +64,15 @@ function OptimalControlProblems.quadrotor(::OptimalControlBackend; nh::Int=60) ## constraints # state constraints - tf ≥ 0.0, (tf_con) - # control constraints + tf ≥ 0.1, (tf_con) -pi / 2 ≤ ϕ(t) ≤ pi / 2, (ϕ_con) -pi / 2 ≤ θ(t) ≤ pi / 2, (θ_con) - cos(θ(t)) * cos(ϕ(t)) ≥ cos(tiltmax), (tiltmax_con) + # control constraints + atmin ≤ at(t) ≤ atmax, (at_con) -dtiltmax ≤ ϕ_dot(t) ≤ dtiltmax, (ϕdot_con) -dtiltmax ≤ θ_dot(t) ≤ dtiltmax, (θdot_con) - atmin ≤ at(t) ≤ atmax, (at_con) + # path constraints + cos(θ(t)) * cos(ϕ(t)) ≥ cos(tiltmax), (tiltmax_con) # initial constraints p1(0) == p0[1], (p1_i) p2(0) == p0[2], (p2_i) @@ -85,33 +97,12 @@ function OptimalControlProblems.quadrotor(::OptimalControlBackend; nh::Int=60) tf + ∫(1e-8 * (ϕ(t)^2 + θ(t)^2 + ψ(t)^2 + at(t)^2) + (1e2 * (ψ(t) - u0[3])^2)) → min end - function dynamics(x, u) - p1, p2, p3, v1, v2, v3, ϕ, θ = x - at, ϕ_dot, θ_dot, ψ = u - - cr = cos(ϕ) - sr = sin(ϕ) - cp = cos(θ) - sp = sin(θ) - cy = cos(ψ) - sy = sin(ψ) - R = [ - (cy*cp) (cy * sp * sr-sy * cr) (cy * sp * cr+sy * sr) - (sy*cp) (sy * sp * sr+cy * cr) (sy * sp * cr-cy * sr) - (-sp) (cp*sr) (cp*cr) - ] - at_ = R * [0; 0; at] - g_ = [0; 0; -g] - a = g_ + at_ - - return [v1, v2, v3, a[1], a[2], a[3], ϕ_dot, θ_dot] - end - # Initial guess - init = (control=[10, 0.0, 0.0, 0.0],) + xinit = t -> [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1] # [p1, p2, p3, v1, v2, v3, ϕ, θ] + uinit = [10.0, 0.1, 0.1, 0.1] # [at, ϕ_dot, θ_dot, ψ] + varinit = [0.1] # [tf] + init = (state=xinit, control=uinit, variable=varinit) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/robbins.jl b/ext/OptimalControlModels/robbins.jl index 199fa744..d8bf2524 100644 --- a/ext/OptimalControlModels/robbins.jl +++ b/ext/OptimalControlModels/robbins.jl @@ -21,10 +21,10 @@ function OptimalControlProblems.robbins(::OptimalControlBackend; nh::Int=100, N: end # Initial guess - init = () + xinit = t -> [0.1, 0.1, 0.1] # [x1, x2, x3] + uinit = [0.1] # [u] + init = (state=xinit, control=uinit) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/robot.jl b/ext/OptimalControlModels/robot.jl index a4000077..83541ac2 100644 --- a/ext/OptimalControlModels/robot.jl +++ b/ext/OptimalControlModels/robot.jl @@ -86,7 +86,5 @@ function OptimalControlProblems.robot(::OptimalControlBackend; nh::Int=100) init = (state=xinit, control=uinit, variable=1.0) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/rocket.jl b/ext/OptimalControlModels/rocket.jl index 4b0e658c..4fd6ddba 100644 --- a/ext/OptimalControlModels/rocket.jl +++ b/ext/OptimalControlModels/rocket.jl @@ -82,7 +82,5 @@ function OptimalControlProblems.rocket(::OptimalControlBackend; nh::Int=100) init = (time=time_vec, state=xinit, control=T_max / 2.0, variable=1) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/space_Shuttle.jl b/ext/OptimalControlModels/space_Shuttle.jl deleted file mode 100644 index 2b6f8a6c..00000000 --- a/ext/OptimalControlModels/space_Shuttle.jl +++ /dev/null @@ -1,195 +0,0 @@ -""" -Space Shuttle Reentry Trajectory Problem: - We want to find the optimal trajectory of a space shuttle reentry. - The objective is to minimize the angle of attack at the terminal point. - The problem is formulated as an OptimalControl model. -""" -function OptimalControlProblems.space_shuttle(::OptimalControlBackend; nh::Int=503) - ## Global variables - w = 203000.0 # weight (lb) - g₀ = 32.174 # acceleration (ft/sec^2) - m = w / g₀ # mass (slug) - - ## Aerodynamic and atmospheric forces on the vehicle - ρ₀ = 0.002378 - hᵣ = 23800.0 - Rₑ = 20902900.0 - μ = 0.14076539e17 - S = 2690.0 - a₀ = -0.20704 - a₁ = 0.029244 - b₀ = 0.07854 - b₁ = -0.61592e-2 - b₂ = 0.621408e-3 - c₀ = 1.0672181 - c₁ = -0.19213774e-1 - c₂ = 0.21286289e-3 - c₃ = -0.10117249e-5 - - ## Initial conditions - h_s = 2.6 # altitude (ft) / 1e5 - ϕ_s = deg2rad(0) # longitude (rad) - θ_s = deg2rad(0) # latitude (rad) - v_s = 2.56 # velocity (ft/sec) / 1e4 - γ_s = deg2rad(-1) # flight path angle (rad) - ψ_s = deg2rad(90) # azimuth (rad) - α_s = deg2rad(0) # angle of attack (rad) - β_s = deg2rad(0) # bank angle (rad) - t_s = 1.00 # time step (sec) - - ## Final conditions, the so-called Terminal Area Energy Management (TAEM) - h_t = 0.8 # altitude (ft) / 1e5 - v_t = 0.25 # velocity (ft/sec) / 1e4 - γ_t = deg2rad(-5) # flight path angle (rad) - tf = 2009.0 # final time (sec) - t0 = 0.0 # initial time (sec) - - ocp = @def begin - ## parameters - w = 203000.0 # weight (lb) - g₀ = 32.174 # acceleration (ft/sec^2) - m = w / g₀ # mass (slug) - ρ₀ = 0.002378 - hᵣ = 23800.0 - Rₑ = 20902900.0 - μ = 0.14076539e17 - S = 2690.0 - a₀ = -0.20704 - a₁ = 0.029244 - b₀ = 0.07854 - b₁ = -0.61592e-2 - b₂ = 0.621408e-3 - c₀ = 1.0672181 - c₁ = -0.19213774e-1 - c₂ = 0.21286289e-3 - c₃ = -0.10117249e-5 - h_s = 2.6 # altitude (ft) / 1e5 - ϕ_s = deg2rad(0) # longitude (rad) - θ_s = deg2rad(0) # latitude (rad) - v_s = 2.56 # velocity (ft/sec) / 1e4 - γ_s = deg2rad(-1) # flight path angle (rad) - ψ_s = deg2rad(90) # azimuth (rad) - α_s = deg2rad(0) # angle of attack (rad) - β_s = deg2rad(0) # bank angle (rad) - t_s = 1.00 - h_t = 0.8 # altitude (ft) / 1e5 - v_t = 0.25 # velocity (ft/sec) / 1e4 - γ_t = deg2rad(-5) # flight path angle (rad) - tf = 2012.0 # final time (sec) - t0 = 0.0 # initial time (sec) - - ## define the problem - t ∈ [t0, tf], time - x ∈ R⁶, state - u ∈ R², control - - ## state variables - scaled_h = x₁ - ϕ = x₂ - θ = x₃ - scaled_v = x₄ - γ = x₅ - ψ = x₆ - - ## control variables - α = u₁ - β = u₂ - - ## constraints - # variable constraints - # state constraints - scaled_h(t) ≥ 0, (scaled_h_con) - deg2rad(-89) ≤ θ(t) ≤ deg2rad(89), (θ_con) - scaled_v(t) ≥ 1e-4, (scaled_v_con) - deg2rad(-89) ≤ γ(t) ≤ deg2rad(89), (γ_con) - # control constraints - deg2rad(-89) ≤ β(t) ≤ deg2rad(1), (β_con) - deg2rad(-90) ≤ α(t) ≤ deg2rad(90), (α_con) - # initial conditions - scaled_h(t0) == h_s, (scaled_h0_con) - ϕ(t0) == ϕ_s, (ϕ0_con) - θ(t0) == θ_s, (θ0_con) - scaled_v(t0) == v_s, (scaled_v0_con) - γ(t0) == γ_s, (γ0_con) - ψ(t0) == ψ_s, (ψ0_con) - # final conditions - scaled_h(tf) == h_t, (scaled_hf_con) - scaled_v(tf) == v_t, (scaled_vf_con) - γ(tf) == γ_t, (γf_con) - - ## dynamics - ẋ(t) == dynamics(x(t), u(t)) - - ## objective - θ(tf) → max - end - - ## dynamics - function dynamics(x, u) - scaled_h, ϕ, θ, scaled_v, γ, ψ = x - α, β = u - h = scaled_h * 1e5 - v = scaled_v * 1e4 - ## Helper functions - c_D = b₀ + b₁ * rad2deg(α) + b₂ * (rad2deg(α)^2) - c_L = a₀ + a₁ * rad2deg(α) - ρ = ρ₀ * exp(-h / hᵣ) - D = (1 / 2) * c_D * S * ρ * (v^2) - L = (1 / 2) * c_L * S * ρ * (v^2) - r = Rₑ + h - g = μ / (r^2) - - ## dynamics - h_dot = v * sin(γ) - ϕ_dot = (v / r) * cos(γ) * sin(ψ) / cos(θ) - θ_dot = (v / r) * cos(γ) * cos(ψ) - v_dot = -(D / m) - g * sin(γ) - γ_dot = (L / (m * v)) * cos(β) + cos(γ) * ((v / r) - (g / v)) - ψ_dot = - (1 / (m * v * cos(γ))) * L * sin(β) + - (v / (r * cos(θ))) * cos(γ) * sin(ψ) * sin(θ) - - return [h_dot, ϕ_dot, θ_dot, v_dot, γ_dot, ψ_dot] - end - - # Initial guess - # Helper function for linear interpolation - function linear_interpolate(x_s, x_t, nh) - return [x_s + (i - 1) / (nh - 1) * (x_t - x_s) for i in 1:nh] - end - # Interpolate each parameter separately - h_interp = linear_interpolate(h_s, h_t, nh) - ϕ_interp = linear_interpolate(ϕ_s, ϕ_s, nh) # no change in longitude - θ_interp = linear_interpolate(θ_s, θ_s, nh) # no change in latitude - v_interp = linear_interpolate(v_s, v_t, nh) - γ_interp = linear_interpolate(γ_s, γ_t, nh) - ψ_interp = linear_interpolate(ψ_s, ψ_s, nh) # no change in azimuth - α_interp = linear_interpolate(α_s, α_s, nh) # no change in angle of attack - β_interp = linear_interpolate(β_s, β_s, nh) # no change in bank angle - t_interp = linear_interpolate(t_s, t_s, nh) # no change in time step - # Combine all interpolated parameters into an array of arrays - interpolated_values = [ - transpose([h, ϕ, θ, v, γ, ψ, α, β, t]) for (h, ϕ, θ, v, γ, ψ, α, β, t) in zip( - h_interp, - ϕ_interp, - θ_interp, - v_interp, - γ_interp, - ψ_interp, - α_interp, - β_interp, - t_interp, - ) - ] - # Create the initial guess by summing the interpolated values - initial_guess = reduce(vcat, interpolated_values) - x_init = [initial_guess[i, 1:6] for i in 1:nh] - u_init = [initial_guess[i, 7:8] for i in 1:nh] - time_vec = LinRange(0.0, t_s * nh * 4, nh) - init = (time=time_vec, state=x_init, control=u_init) - - # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp -end diff --git a/ext/OptimalControlModels/space_shuttle.jl b/ext/OptimalControlModels/space_shuttle.jl new file mode 100644 index 00000000..df1ffb37 --- /dev/null +++ b/ext/OptimalControlModels/space_shuttle.jl @@ -0,0 +1,133 @@ +""" +Space Shuttle Reentry Trajectory Problem: + We want to find the optimal trajectory of a space shuttle reentry. + The objective is to maximize the latitude (cross range) at the terminal point. + The original problem formulated as a JuMP model can be found [here](https://jump.dev/JuMP.jl/stable/tutorials/nonlinear/space_shuttle_reentry_trajectory/) + Note: no heating limit path constraint +""" +function OptimalControlProblems.space_shuttle(::OptimalControlBackend; nh::Int=503) + ## Global variables + w = 203000.0 # weight (lb) + g₀ = 32.174 # acceleration (ft/sec^2) + m = w / g₀ # mass (slug) + + ## Aerodynamic and atmospheric forces on the vehicle + ρ₀ = 0.002378 + hᵣ = 23800.0 + Rₑ = 20902900.0 + μ = 0.14076539e17 + S = 2690.0 + a₀ = -0.20704 + a₁ = 0.029244 + b₀ = 0.07854 + b₁ = -0.61592e-2 + b₂ = 0.621408e-3 + + ## Initial conditions + h_s = 2.6 # altitude (ft) / 1e5 + ϕ_s = deg2rad(0) # longitude (rad) + θ_s = deg2rad(0) # latitude (rad) + v_s = 2.56 # velocity (ft/sec) / 1e4 + γ_s = deg2rad(-1) # flight path angle (rad) + ψ_s = deg2rad(90) # azimuth (rad) + α_s = deg2rad(0) # angle of attack (rad) + β_s = deg2rad(0) # bank angle (rad) + + ## Final conditions, the so-called Terminal Area Energy Management (TAEM) + h_t = 0.8 # altitude (ft) / 1e5 + v_t = 0.25 # velocity (ft/sec) / 1e4 + γ_t = deg2rad(-5) # flight path angle (rad) + + ## dynamics + function dynamics(x, u) + scaled_h, ϕ, θ, scaled_v, γ, ψ = x + α, β = u + h = scaled_h * 1e5 + v = scaled_v * 1e4 + ## Helper functions + c_D = b₀ + b₁ * rad2deg(α) + b₂ * (rad2deg(α)^2) + c_L = a₀ + a₁ * rad2deg(α) + ρ = ρ₀ * exp(-h / hᵣ) + D = (1 / 2) * c_D * S * ρ * (v^2) + L = (1 / 2) * c_L * S * ρ * (v^2) + r = Rₑ + h + g = μ / (r^2) + + ## dynamics + h_dot = v * sin(γ) + ϕ_dot = (v / r) * cos(γ) * sin(ψ) / cos(θ) + θ_dot = (v / r) * cos(γ) * cos(ψ) + v_dot = -(D / m) - g * sin(γ) + γ_dot = (L / (m * v)) * cos(β) + cos(γ) * ((v / r) - (g / v)) + ψ_dot = + (1 / (m * v * cos(γ))) * L * sin(β) + + (v / (r * cos(θ))) * cos(γ) * sin(ψ) * sin(θ) + + return [h_dot / 1e5, ϕ_dot, θ_dot, v_dot / 1e4, γ_dot, ψ_dot] + end + + ocp = @def begin + + ## define the problem + tf ∈ R¹, variable + t ∈ [0, tf], time + x ∈ R⁶, state + u ∈ R², control + + ## state variables + scaled_h = x₁ + ϕ = x₂ + θ = x₃ + scaled_v = x₄ + γ = x₅ + ψ = x₆ + + ## control variables + α = u₁ + β = u₂ + + ## constraints + 1750 ≤ tf ≤ 2250 # NB jump with 503 steps between 3.5 and 4.5 + # state constraints + 0 ≤ scaled_h(t) ≤ Inf, (scaled_h_con) + deg2rad(-89) ≤ θ(t) ≤ deg2rad(89), (θ_con) + 0 ≤ scaled_v(t) ≤ Inf, (scaled_v_con) + deg2rad(-89) ≤ γ(t) ≤ deg2rad(89), (γ_con) + # control constraints + deg2rad(-90) ≤ α(t) ≤ deg2rad(90), (α_con) + deg2rad(-89) ≤ β(t) ≤ deg2rad(1), (β_con) + + # initial conditions + scaled_h(0) == h_s, (scaled_h0_con) + ϕ(0) == ϕ_s, (ϕ0_con) + θ(0) == θ_s, (θ0_con) + scaled_v(0) == v_s, (scaled_v0_con) + γ(0) == γ_s, (γ0_con) + ψ(0) == ψ_s, (ψ0_con) + # final conditions + scaled_h(tf) == h_t, (scaled_hf_con) + scaled_v(tf) == v_t, (scaled_vf_con) + γ(tf) == γ_t, (γf_con) + + ## dynamics + ẋ(t) == dynamics(x(t), u(t)) + + ## objective + θ(tf) → max + end + + # initial guess: linear interpolation for h, v, gamma (NB. t0 = 0), constant for the rest + # variable time step seems to be initialized at 1 in jump + # note that ipopt will project the initial guess inside the bounds anyway. + tf_init = 500 + x_init = t -> [ h_s + t / tf_init * (h_t - h_s) , + ϕ_s, + θ_s, + v_s + t / tf_init * (v_t - v_s), + γ_s + t / tf_init * (γ_t - γ_s), + ψ_s] + init = (state=x_init, control=[α_s, β_s], variable=[tf_init]) + + # NLPModel + DOCP + return direct_transcription(ocp; init=init, grid_size=nh) +end diff --git a/ext/OptimalControlModels/steering.jl b/ext/OptimalControlModels/steering.jl index 976c1259..5639902b 100644 --- a/ext/OptimalControlModels/steering.jl +++ b/ext/OptimalControlModels/steering.jl @@ -63,7 +63,5 @@ function OptimalControlProblems.steering(::OptimalControlBackend; nh::Int=100) init = (state=xinit, control=0.0, variable=1.0) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/truck_trailer.jl b/ext/OptimalControlModels/truck_trailer.jl index 0c158661..d9e462fd 100644 --- a/ext/OptimalControlModels/truck_trailer.jl +++ b/ext/OptimalControlModels/truck_trailer.jl @@ -7,12 +7,9 @@ The Truck Trailer Problem: function OptimalControlProblems.truck_trailer( ::OptimalControlBackend; data::Array{Float64,2}=[0.4 0.1 0.2; 1.1 0.2 0.2; 0.8 0.1 0.2], - nh::Int=100, + nh::Int=250, ) # parameters - if size(data) != (3, 3) - error("The input matrix must be 3x3.") - end L0 = data[1, 1] M0 = data[1, 2] W0 = data[1, 3] @@ -30,32 +27,28 @@ function OptimalControlProblems.truck_trailer( theta0_t0 = 0.0 x2_tf = 0.0 y2_tf = -2.0 - theta2_tf = 2 * pi / 4 - theta1_tf = 2 * pi / 4 - theta0_tf = 2 * pi / 4 + theta2_tf = pi / 2 + theta1_tf = pi / 2 + theta0_tf = pi / 2 + + function dynamics(x, u) + x2, y2, theta0, theta1, theta2, v0, delta0 = x + d_v0, d_delta0 = u + + beta01 = theta0 - theta1 + beta12 = theta1 - theta2 + dtheta0 = v0 / L0 * tan(delta0) + dtheta1 = v0 / L1 * sin(beta01) - M0 / L1 * cos(beta01) * dtheta0 + v1 = v0 * cos(beta01) + M0 * sin(beta01) * dtheta0 + dtheta2 = v1 / L2 * sin(beta12) - M1 / L2 * cos(beta12) * dtheta1 + v2 = v1 * cos(beta12) + M1 * sin(beta12) * dtheta1 + dx2 = v2 * cos(theta2) + dy2 = v2 * sin(theta2) + + return [dx2, dy2, dtheta0, dtheta1, dtheta2, d_v0, d_delta0] + end ocp = @def begin - # parameters - L0 = data[1, 1] - M0 = data[1, 2] - W0 = data[1, 3] - L1 = data[2, 1] - M1 = data[2, 2] - W1 = data[2, 3] - L2 = data[3, 1] - M2 = data[3, 2] - W2 = data[3, 3] - speedf = 1 - x2_t0 = 0.0 - y2_t0 = 0.0 - theta2_t0 = 0.0 - theta1_t0 = 0.0 - theta0_t0 = 0.0 - x2_tf = 0.0 - y2_tf = -2.0 - theta2_tf = 2 * pi / 4 - theta1_tf = 2 * pi / 4 - theta0_tf = 2 * pi / 4 ## define the problem tf ∈ R, variable @@ -83,16 +76,22 @@ function OptimalControlProblems.truck_trailer( y0 = y_1 + L1 * sin(theta1) + M0 * sin(theta0) ## constraints + 1.0 ≤ tf ≤ 1000 + # state constraints -pi / 2 ≤ theta0(t) ≤ pi / 2, (theta0_con) -pi / 2 ≤ theta1(t) ≤ pi / 2, (theta1_con) + # control constraints -0.2 * speedf ≤ v0(t) ≤ 0.2 * speedf, (v0_con) -pi / 6 ≤ delta0(t) ≤ pi / 6, (delta0_con) + + # path constraints -1 ≤ d_v0(t) ≤ 1, (v0_dot_con) -pi / 10 ≤ d_delta0(t) ≤ pi / 10, (delta0_dot_con) -pi / 2 ≤ beta01(t) ≤ pi / 2, (beta01_con) -pi / 2 ≤ beta12(t) ≤ pi / 2, (beta12_con) + # initial conditions x_2(0) == x2_t0, (x2_t0_con) y_2(0) == y2_t0, (y2_t0_con) @@ -113,28 +112,12 @@ function OptimalControlProblems.truck_trailer( tf + ∫((beta01(t)^2) + (beta12(t)^2)) → min end - function dynamics(x, u) - x2, y2, theta0, theta1, theta2, v0, delta0 = x - d_v0, d_delta0 = u - - beta01 = theta0 - theta1 - beta12 = theta1 - theta2 - dtheta0 = v0 / L0 * tan(delta0) - dtheta1 = v0 / L1 * sin(beta01) - M0 / L1 * cos(beta01) * dtheta0 - v1 = v0 * cos(beta01) + M0 * sin(beta01) * dtheta0 - dtheta2 = v1 / L2 * sin(beta12) - M1 / L2 * cos(beta12) * dtheta1 - v2 = v1 * cos(beta12) + M1 * sin(beta12) * dtheta1 - dx2 = v2 * cos(theta2) - dy2 = v2 * sin(theta2) - - return [dx2, dy2, dtheta0, dtheta1, dtheta2, d_v0, d_delta0] - end - # Initial guess - init = (state=[0, 0, 0.1, 0.0, 0.0, -0.2, 0],) + xinit = t -> [0.1, 0.1, 0.1, 0.0, 0.0, -0.2, 0.1] # [x2, y2, theta0, theta1, theta2, v0, delta0] + uinit = [0.1, 0.1] # [d_v0, d_delta0] + varinit = [10] # [tf] + init = (state=xinit, control=uinit, variable=varinit) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/ext/OptimalControlModels/vanderpol.jl b/ext/OptimalControlModels/vanderpol.jl index 250e8fce..58501d48 100644 --- a/ext/OptimalControlModels/vanderpol.jl +++ b/ext/OptimalControlModels/vanderpol.jl @@ -18,10 +18,10 @@ function OptimalControlProblems.vanderpol(::OptimalControlBackend; nh::Int=100) end # Initial guess - init = () + xinit = t -> [0.1, 0.1] # [x1, x2] + uinit = [0.1] # [u] + init = (state=xinit, control=uinit) # NLPModel + DOCP - docp, nlp = direct_transcription(ocp; init=init, grid_size=nh) - - return docp, nlp + return direct_transcription(ocp; init=init, grid_size=nh) end diff --git a/src/OptimalControlProblems.jl b/src/OptimalControlProblems.jl index 6bd27108..f4f1e75b 100644 --- a/src/OptimalControlProblems.jl +++ b/src/OptimalControlProblems.jl @@ -37,6 +37,10 @@ const infos = [ :nvar :ncon :minimize + :state_name + :costate_name + :control_name + :time ] const types = [ @@ -45,6 +49,10 @@ const types = [ Union{Int,Nothing}, Union{Int,Nothing}, Union{Bool,Nothing}, + Union{Vector{String}, String, Nothing}, + Union{Vector{String}, String, Nothing}, + Union{Vector{String}, String, Nothing}, + Union{Tuple{String, String, Union{Real,Nothing}}, Nothing}, ] """ @@ -71,6 +79,33 @@ for i in 1:number_of_problems end end -export JuMPBackend, OptimalControlBackend +# ------- Available Problems Function ------- +""" + available_problems() + +Returns the list of problems that are currently working based on the last test execution. +The list is read from a cache file that gets updated every time tests are run. +If no cache exists, returns an empty list with a warning to run tests first. +""" +function available_problems() + cache_file = joinpath(dirname(@__FILE__), "..", "available_problems_cache.txt") + if isfile(cache_file) + try + content = read(cache_file, String) + if !isempty(strip(content)) + # Parse symbols from the file + lines = split(strip(content), '\n') + return [Symbol(strip(line)) for line in lines if !isempty(strip(line))] + end + catch e + @warn "Error reading the cache: $e" + end + end + # Default list if the cache does not exist or is empty + @warn "Available problems cache not found. Run the tests to update the list." + return Symbol[] +end + +export JuMPBackend, OptimalControlBackend, available_problems end diff --git a/test/Project.toml b/test/Project.toml index cc8a3803..ad774274 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,16 +1,22 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71" OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.8" CTBase = "0.16" +Interpolations = "0.15, 0.16" Ipopt = "1.9" JuMP = "1.25" NLPModelsIpopt = "0.10" OptimalControl = "1.0" +Plots = "1.40" +Test = "1" diff --git a/test/figures/beam.png b/test/figures/beam.png new file mode 100644 index 00000000..789f5e4f Binary files /dev/null and b/test/figures/beam.png differ diff --git a/test/figures/bioreactor.png b/test/figures/bioreactor.png new file mode 100644 index 00000000..86a199c6 Binary files /dev/null and b/test/figures/bioreactor.png differ diff --git a/test/figures/cart_pendulum.png b/test/figures/cart_pendulum.png new file mode 100644 index 00000000..fcac1b1c Binary files /dev/null and b/test/figures/cart_pendulum.png differ diff --git a/test/figures/chain.png b/test/figures/chain.png new file mode 100644 index 00000000..e57d2554 Binary files /dev/null and b/test/figures/chain.png differ diff --git a/test/figures/dielectrophoretic_particle.png b/test/figures/dielectrophoretic_particle.png new file mode 100644 index 00000000..449cd7b5 Binary files /dev/null and b/test/figures/dielectrophoretic_particle.png differ diff --git a/test/figures/double_oscillator.png b/test/figures/double_oscillator.png new file mode 100644 index 00000000..6e0b0240 Binary files /dev/null and b/test/figures/double_oscillator.png differ diff --git a/test/figures/ducted_fan.png b/test/figures/ducted_fan.png new file mode 100644 index 00000000..97d7ea48 Binary files /dev/null and b/test/figures/ducted_fan.png differ diff --git a/test/figures/electric_vehicle.png b/test/figures/electric_vehicle.png new file mode 100644 index 00000000..a782f349 Binary files /dev/null and b/test/figures/electric_vehicle.png differ diff --git a/test/figures/glider.png b/test/figures/glider.png new file mode 100644 index 00000000..00f9c7f3 Binary files /dev/null and b/test/figures/glider.png differ diff --git a/test/figures/insurance.png b/test/figures/insurance.png new file mode 100644 index 00000000..de970ee0 Binary files /dev/null and b/test/figures/insurance.png differ diff --git a/test/figures/jackson.png b/test/figures/jackson.png new file mode 100644 index 00000000..92e6d079 Binary files /dev/null and b/test/figures/jackson.png differ diff --git a/test/figures/moonlander.png b/test/figures/moonlander.png new file mode 100644 index 00000000..2f72c330 Binary files /dev/null and b/test/figures/moonlander.png differ diff --git a/test/figures/quadrotor.png b/test/figures/quadrotor.png new file mode 100644 index 00000000..9b0160ea Binary files /dev/null and b/test/figures/quadrotor.png differ diff --git a/test/figures/robbins.png b/test/figures/robbins.png new file mode 100644 index 00000000..adfedabc Binary files /dev/null and b/test/figures/robbins.png differ diff --git a/test/figures/robot.png b/test/figures/robot.png new file mode 100644 index 00000000..708ea7cf Binary files /dev/null and b/test/figures/robot.png differ diff --git a/test/figures/rocket.png b/test/figures/rocket.png new file mode 100644 index 00000000..f1f54b20 Binary files /dev/null and b/test/figures/rocket.png differ diff --git a/test/figures/steering.png b/test/figures/steering.png new file mode 100644 index 00000000..2fd70ae0 Binary files /dev/null and b/test/figures/steering.png differ diff --git a/test/figures/truck_trailer.png b/test/figures/truck_trailer.png new file mode 100644 index 00000000..1bedac26 Binary files /dev/null and b/test/figures/truck_trailer.png differ diff --git a/test/figures/vanderpol.png b/test/figures/vanderpol.png new file mode 100644 index 00000000..bdc89c50 Binary files /dev/null and b/test/figures/vanderpol.png differ diff --git a/test/runtests.jl b/test/runtests.jl index e53a05ec..3f5a15ca 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,27 +6,86 @@ using NLPModelsIpopt using OptimalControl using OptimalControlProblems using Test +using Plots +using Interpolations +include("utils.jl") -# Parameters -const tol = 1e-6 +# Parameters for the solvers +const tol = 1e-8 const mu_strategy = "adaptive" const sb = "yes" -const constr_viol_tol = 1e-6 const max_iter = 1000 const max_wall_time = 500.0 -# -@testset verbose = true showtiming = true "OptimalControlProblems tests" begin +# Collect all the problem from OptimalControlProblems +all_names = names(OptimalControlProblems; all=true) +functions_list = filter( + x -> + isdefined(OptimalControlProblems, x) && + isa(getfield(OptimalControlProblems, x), Function) && + !startswith(string(x), "#") && + !(x in [:eval, :include, :available_problems]), + all_names, +) + +# Remove from the tests the problems that are not working +pbs_with_issues = [ +] +functions_list = setdiff(functions_list, pbs_with_issues) + +# The list of all the problems to test +const list_of_problems = deepcopy(functions_list) + +# The final list of problems for which the tests pass +list_of_problems_final = deepcopy(functions_list) + +# Tests +const verbose = true # print or not details during tests +@testset "OptimalControlProblems tests" verbose=verbose showtiming=true begin + for name in ( - :aqua, - :JuMP, - :OptimalControl, + #:aqua, + #:JuMP, # convergence tests for JuMP models + #:OptimalControl, # convergence tests for OptimalControl models + :Comparison, # comparison between OptimalControl and JuMP + #:quick, # quick comparison: objective rel error only ) - @testset "$(name)" begin + @testset "$(name)" verbose=verbose begin test_name = Symbol(:test_, name) println("Testing: " * string(name)) include("$(test_name).jl") @eval $test_name() end end + + # note: not sure I undertand the purpose of the cache ? + # the list of available problems is not up to date when the tests are run + # and we need to run the tests another time to get the proper list + @testset "available_problems" verbose=verbose begin + if list_of_problems_final == available_problems() + @test list_of_problems_final == available_problems() + else + println("Available problems: ", available_problems()) + println("Passed problems: ",list_of_problems_final) + @test list_of_problems_final == available_problems() broken=true + end + end + + +end + +# +println("List of problems working:", list_of_problems_final) + +# Sauvegarder la liste des problèmes fonctionnels dans un fichier de cache +cache_file = joinpath(@__DIR__, "..", "available_problems_cache.txt") +try + open(cache_file, "w") do f + for problem in list_of_problems_final + println(f, string(problem)) + end + end + println("Cache des problèmes disponibles mis à jour: $cache_file") +catch e + @warn "Impossible de sauvegarder le cache des problèmes: $e" end diff --git a/test/test_Comparison.jl b/test/test_Comparison.jl new file mode 100644 index 00000000..261faf56 --- /dev/null +++ b/test/test_Comparison.jl @@ -0,0 +1,340 @@ +function test_Comparison() + + test_init = false + + # Comparison Parameters + ε = 1e-2 + p = 2 + + # options for solvers + kwargs_init = Dict( + :print_level => 0, + :tol => tol, + :mu_strategy => mu_strategy, + :sb => sb, + :max_iter => 0, + :max_wall_time => max_wall_time, + ) + + kwargs = Dict( + :print_level => 0, + :tol => tol, + :mu_strategy => mu_strategy, + :sb => sb, + :max_iter => max_iter, + :max_wall_time => max_wall_time, + ) + + println() + println("\033[1m###########################\033[0m") + println("\033[1m####### COMPARISON ########\033[0m") + println("\033[1m###########################\033[0m") + println() + + for f in list_of_problems + nh = OptimalControlProblems.metadata[f][:nh] + @testset "$(string(f))" verbose = verbose begin + println("############ TEST $f #############") + println() + + if test_init + #================== INIT ======================# + + ########## OptimalControl ########## + + # Set up the OptimalControl model + docp_init, OC_model_init = OptimalControlProblems.eval(f)(OptimalControlBackend()) # +++ UPDATE + + # Solve the problem + nlp_sol_init = NLPModelsIpopt.ipopt(OC_model_init; kwargs_init...) + + # Build the solution + sol_init = build_OCP_solution(docp_init; primal=nlp_sol_init.solution, dual=nlp_sol_init.multipliers) + + # Retrieves values of variables + x_init_oc = state(sol_init) + p_init_oc = costate(sol_init) + u_init_oc = control(sol_init) + + ############### JuMP ############### + + # Set up the JuMP model + JuMP_init_model = OptimalControlProblems.eval(f)(JuMPBackend()) + set_optimizer(JuMP_init_model, Ipopt.Optimizer) + set_silent(JuMP_init_model) + set_optimizer_attribute(JuMP_init_model, "tol", tol) + set_optimizer_attribute(JuMP_init_model, "max_iter", 0) + set_optimizer_attribute(JuMP_init_model, "mu_strategy", mu_strategy) + set_optimizer_attribute(JuMP_init_model, "linear_solver", "mumps") + set_optimizer_attribute(JuMP_init_model, "max_wall_time", max_wall_time) + set_optimizer_attribute(JuMP_init_model, "sb", sb) + + # Solve the model + optimize!(JuMP_init_model) + + # Retrieves values of variables + time_data, time_var_name, time_value = OptimalControlProblems.metadata[f][:time] + if time_data == "final_time" + if time_value !== nothing + tf_init = time_value + else + tf_init = value.(JuMP_init_model[Symbol(time_var_name)]) + end + h_init = tf_init / nh + elseif time_data == "step" + if time_value !== nothing + h_init = time_value + else + h_init = value.(JuMP_init_model[Symbol(time_var_name)]) + end + tf_init = h_init * nh + end + t_init = Vector((0:nh) * h_init) + + x_vars = OptimalControlProblems.metadata[f][:state_name] + p_vars = OptimalControlProblems.metadata[f][:costate_name] + u_vars = OptimalControlProblems.metadata[f][:control_name] + + x_jmp_vars_init = [JuMP.value.(JuMP_init_model[Symbol(xv)]) for xv in x_vars] + inds_x_init = axes(x_jmp_vars_init[1], 1) + x_jmp_init = [[x_jmp_vars_init[j][i] for j in 1:length(x_vars)] for i in inds_x_init] + + u_jmp_vars_init = [JuMP.value.(JuMP_init_model[Symbol(uv)]) for uv in u_vars] + inds_u_init = axes(u_jmp_vars_init[1], 1) + u_jmp_init = [[u_jmp_vars_init[j][i] for j in 1:length(u_vars)] for i in inds_u_init] + + p_jmp_vars_init = [JuMP.dual.(JuMP_init_model[Symbol(pv)]) for pv in p_vars] + inds_p_init = axes(p_jmp_vars_init[1], 1) + p_jmp_init = -[[p_jmp_vars_init[j][i] for j in 1:length(p_vars)] for i in inds_p_init] + p_jmp_init = costateInterpolation(p_jmp_init, t_init) + + ############ TEST ############ + @testset "init" verbose=verbose begin + print("Init:\n") + for k in 1:length(x_jmp_init[1]) + dist_x_init = abs(x_init_oc(0)[k] - x_jmp_init[1][k]) + print(" Test x$k : ") + @testset "x$k" verbose=verbose begin + if !(dist_x_init < ε) + print("$dist_x_init < $ε \033[1;33mTest Broken\033[0m\n") + @test dist_x_init < ε broken=true + global list_of_problems_final + list_of_problems_final = setdiff(list_of_problems_final, [f]) + else + print("$dist_x_init < $ε \033[1;32mTest Passed\033[0m\n") + @test dist_x_init < ε + end + end + end + + for k in 1:length(p_jmp_init[1]) + dist_p_init = abs(p_init_oc(0)[k] - p_jmp_init[1][k]) + print(" Test p$k : ") + @testset "p$k" verbose=verbose begin + if !(dist_p_init < ε) + print("$dist_p_init < $ε \033[1;33mTest Broken\033[0m\n") + @test dist_p_init < ε broken=true + else + print("$dist_p_init < $ε \033[1;32mTest Passed\033[0m\n") + @test dist_p_init < ε + end + end + end + + for k in 1:length(u_jmp_init[1]) + dist_u_init = abs(u_init_oc(0)[k] - u_jmp_init[1][k]) + print(" Test u$k : ") + @testset "u$k" verbose=verbose begin + if !(dist_u_init < ε) + print("$dist_u_init < $ε \033[1;33mTest Broken\033[0m\n") + @test dist_u_init < ε broken=true + global list_of_problems_final + list_of_problems_final = setdiff(list_of_problems_final, [f]) + else + print("$dist_u_init < $ε \033[1;32mTest Passed\033[0m\n") + @test dist_u_init < ε + end + end + end + end + end + + #======================= END INIT ========================# + + #======================= Lp + OBJECTIVE =================# + + ########## OptimalControl ########## + + # Set up the OptimalControl model + docp, OC_model = OptimalControlProblems.eval(f)(OptimalControlBackend()) + + # Solve the problem + nlp_sol = NLPModelsIpopt.ipopt(OC_model; kwargs...) + + # Build the solution + sol = build_OCP_solution(docp; primal=nlp_sol.solution, dual=nlp_sol.multipliers) + + # Retrieves values of variables + x_oc = state(sol) + p_oc = costate(sol) + u_oc = control(sol) + obj_oc = objective(sol) + + ############### JuMP ############### + + # Set up the JuMP model + JuMP_model = OptimalControlProblems.eval(f)(JuMPBackend()) + set_optimizer(JuMP_model, Ipopt.Optimizer) + set_silent(JuMP_model) + set_optimizer_attribute(JuMP_model, "tol", tol) + set_optimizer_attribute(JuMP_model, "max_iter", max_iter) + set_optimizer_attribute(JuMP_model, "mu_strategy", mu_strategy) + set_optimizer_attribute(JuMP_model, "linear_solver", "mumps") + set_optimizer_attribute(JuMP_model, "max_wall_time", max_wall_time) + set_optimizer_attribute(JuMP_model, "sb", sb) + + # Solve the model + optimize!(JuMP_model) + + time_data, time_var_name, time_value = OptimalControlProblems.metadata[f][:time] + if time_data == "final_time" + if time_value !== nothing + tf = time_value + else + tf = value.(JuMP_model[Symbol(time_var_name)]) + end + h = tf / nh + elseif time_data == "step" + if time_value !== nothing + h = time_value + else + h = value.(JuMP_model[Symbol(time_var_name)]) + end + tf = h * nh + end + t = Vector((0:nh) * h) + + x_jmp_vars = [JuMP.value.(JuMP_model[Symbol(xv)]) for xv in x_vars] + inds_x = axes(x_jmp_vars[1], 1) + x_jmp = [[x_jmp_vars[j][i] for j in 1:length(x_vars)] for i in inds_x] + + u_jmp_vars = [JuMP.value.(JuMP_model[Symbol(uv)]) for uv in u_vars] + inds_u = axes(u_jmp_vars[1], 1) + u_jmp = [[u_jmp_vars[j][i] for j in 1:length(u_vars)] for i in inds_u] + + p_jmp_vars = [JuMP.dual.(JuMP_model[Symbol(pv)]) for pv in p_vars] + inds_p = axes(p_jmp_vars[1], 1) + p_jmp = -[[p_jmp_vars[j][i] for j in 1:length(p_vars)] for i in inds_p] + p_jmp = costateInterpolation(p_jmp, t) + + obj_jmp = objective_value(JuMP_model) + + dist_obj = abs(obj_oc - obj_jmp) / (obj_oc + obj_jmp) / 2. + @testset "objective" verbose=verbose begin + print("Objective:\n") + print(" Test objective : ") + if !(dist_obj < ε) + print("$dist_obj < $ε \033[1;33mTest Broken\033[0m\n") + println("Jump ", obj_jmp, " vs OC ", obj_oc) + @test dist_obj < ε broken=true + global list_of_problems_final + list_of_problems_final = setdiff(list_of_problems_final, [f]) + else + print("$dist_obj < $ε \033[1;32mTest Passed\033[0m\n") + @test dist_obj < ε + end + end + + plots_x = Vector{Any}() + plots_p = Vector{Any}() + plots_u = Vector{Any}() + + # +++ use relative error + @testset "norm_L$p" verbose=verbose begin + print("Norm_L$p:\n") + for k in 1:length(x_jmp[1]) + dist_x_Lp = norm_Lp([x_oc((i - 1) * h)[k] - x_jmp[i][k] for i in 1:nh+1], p, h) + print(" Test x$k : ") + @testset "x$k" verbose = verbose begin + if !(dist_x_Lp < ε) + print("$dist_x_Lp < $ε \033[1;33mTest Broken\033[0m\n") + @test dist_x_Lp < ε broken=true + global list_of_problems_final + list_of_problems_final = setdiff(list_of_problems_final, [f]) + else + print("$dist_x_Lp < $ε \033[1;32mTest Passed\033[0m\n") + @test dist_x_Lp < ε + end + end + px = plot(plot(sol)[k]; line=2, label="OptimalControl") # OptimalControl + px = plot!(t, [x_jmp[i][k] for i in 1:nh+1]; xlabel="t", ylabel=string(x_vars[k]), legend=false, line=2, color="red", linestyle=:dash, label="JuMP") # JuMP + push!(plots_x, px) + end + + for k in 1:length(p_jmp[1]) + dist_p_Lp = norm_Lp([p_oc((i - 1) * h)[k] - p_jmp[i][k] for i in 1:nh+1], p, h) + print(" Test p$k : ") + @testset "p$k" verbose = verbose begin + if !(dist_p_Lp < ε) + print("$dist_p_Lp < $ε \033[1;33mTest Broken\033[0m\n") + @test dist_p_Lp < ε broken=true + global list_of_problems_final + list_of_problems_final = setdiff(list_of_problems_final, [f]) + else + print("$dist_p_Lp < $ε \033[1;32mTest Passed\033[0m\n") + @test dist_p_Lp < ε + end + end + pp = plot(plot(sol)[length(x_jmp[1])+k]; line=2, label="OptimalControl") # OptimalControl + pp = plot!(t, [p_jmp[i][k] for i in 1:nh+1]; xlabel="t", ylabel="p_" * string(x_vars[k]), legend=false, line=2, color="red", linestyle=:dash, label="JuMP") # JuMP + push!(plots_p, pp) + end + + for k in 1:length(u_jmp[1]) + dist_u_Lp = norm_Lp([u_oc((i - 1) * h)[k] - u_jmp[i][k] for i in 1:nh+1], p, h) + print(" Test u$k : ") + @testset "u$k" verbose = verbose begin + if !(dist_u_Lp < ε) + print("$dist_u_Lp < $ε \033[1;33mTest Broken\033[0m\n") + @test dist_u_Lp < ε broken=true + global list_of_problems_final + list_of_problems_final = setdiff(list_of_problems_final, [f]) + else + print("$dist_u_Lp < $ε \033[1;32mTest Passed\033[0m\n") + @test dist_u_Lp < ε + end + end + pu = plot(plot(sol)[length(x_jmp[1])+length(p_jmp[1])+k]; line=2, label="OptimalControl") # OptimalControl + pu = plot!(t, [u_jmp[i][k] for i in 1:nh+1]; xlabel="t", ylabel=string(u_vars[k]), legend=false, line=2, color="red", linestyle=:dash, label="JuMP") # JuMP + push!(plots_u, pu) + end + end + + all_plots = vcat(plots_x, plots_p, plots_u) + n = length(all_plots) + + figdir = joinpath(@__DIR__, "figures") + isdir(figdir) || mkpath(figdir) + + fig = plot(all_plots..., layout=(n, 1), size=(900, 220 * n), suptitle="Comparison $f", titlefont=font(18)) + display(fig) + + savefig(fig, joinpath(figdir, "$f" * ".png")) + + println() + println("############ END TEST $f #############") + println() + + #================== END Lp + OBJECTIVE =====================# + + end + end + + println() + println("\033[1m###########################\033[0m") + println("\033[1m##### END COMPARISON ######\033[0m") + println("\033[1m###########################\033[0m") + println() + + +end \ No newline at end of file diff --git a/test/test_JuMP.jl b/test/test_JuMP.jl index 7ed6f691..259af4e8 100644 --- a/test/test_JuMP.jl +++ b/test/test_JuMP.jl @@ -1,28 +1,20 @@ # test_JuMP_optimality function test_JuMP() - # Collecting all the OptimalControlProblems.JuMPModels models - all_names = names(OptimalControlProblems; all=true) - functions_list = filter( - x -> - isdefined(OptimalControlProblems, x) && - isa(getfield(OptimalControlProblems, x), Function) && - !startswith(string(x), "#") && - !(x in [:eval, :include]), - all_names, - ) - pbs_with_issues = [:cart_pendulum, :truck_trailer] - functions_list = setdiff(functions_list, pbs_with_issues) + println() + println("\033[1m###########################\033[0m") + println("\033[1m### TEST CONVERGED JuMP ###\033[0m") + println("\033[1m###########################\033[0m") + println() - for f in functions_list - @testset "$(f)" begin + for f in list_of_problems + @testset "$(f)" verbose=verbose begin println(" $f:") # Set up the model model = OptimalControlProblems.eval(f)(JuMPBackend()) set_optimizer(model, Ipopt.Optimizer) set_silent(model) set_optimizer_attribute(model, "tol", tol) - set_optimizer_attribute(model, "constr_viol_tol", constr_viol_tol) set_optimizer_attribute(model, "max_iter", max_iter) set_optimizer_attribute(model, "mu_strategy", mu_strategy) set_optimizer_attribute(model, "linear_solver", "mumps") @@ -32,12 +24,23 @@ function test_JuMP() print(" First solve: "); @time optimize!(model) print(" Second solve: "); @time optimize!(model) # Test that the solver found an optimal solution - println(" termination_status = $(termination_status(model))\n") - if f == :truck_trailer || f == :quadrotor - @test (termination_status(model) == MOI.LOCALLY_INFEASIBLE) || (termination_status(model) == MOI.ITERATION_LIMIT) - else + println(" termination_status = $(termination_status(model)) objective = $(objective_value(model))\n") + if termination_status(model) == MOI.LOCALLY_SOLVED @test termination_status(model) == MOI.LOCALLY_SOLVED + print("JuMP: $f converged : \033[1;32mTest Passed\033[0m\n") + else + @test termination_status(model) == MOI.LOCALLY_SOLVED broken=true + print("JuMP : $f converged : \033[1;33mTest Broken\033[0m\n") + global list_of_problems_final + list_of_problems_final = setdiff(list_of_problems_final, [f]) end end end + + println() + println("\033[1m###############################\033[0m") + println("\033[1m### END TEST CONVERGED JuMP ###\033[0m") + println("\033[1m###############################\033[0m") + println() + end diff --git a/test/test_OptimalControl.jl b/test/test_OptimalControl.jl index 698719ec..848c6443 100644 --- a/test/test_OptimalControl.jl +++ b/test/test_OptimalControl.jl @@ -1,44 +1,48 @@ # test_OptimalControl_optimality function test_OptimalControl() - # Collecting all the OptimalControlProblems.OptimalControlModels models - all_names = names(OptimalControlProblems; all=true) - functions_list = filter( - x -> - isdefined(OptimalControlProblems, x) && - isa(getfield(OptimalControlProblems, x), Function) && - !startswith(string(x), "#") && - !(x in [:eval, :include]), - all_names, - ) - - pbs_with_issues = [:glider, :moonlander] - functions_list = setdiff(functions_list, pbs_with_issues) kwargs = Dict( :print_level => 0, :tol => tol, :mu_strategy => mu_strategy, :sb => sb, - :constr_viol_tol => constr_viol_tol, :max_iter => max_iter, :max_wall_time => max_wall_time, ) - for f in functions_list + println() + println("\033[1m#########################################\033[0m") + println("\033[1m##### TEST CONVERGED OptimalControl #####\033[0m") + println("\033[1m#########################################\033[0m") + println() + + for f in list_of_problems println(" $f:") - @testset "$(f)" begin + @testset "$(f)" verbose=verbose begin # Set up the model - _, model = OptimalControlProblems.eval(f)(OptimalControlBackend()) + _, model = OptimalControlProblems.eval(f)(OptimalControlBackend()) # !+++ UPDATE print(" First solve: "); @time sol = NLPModelsIpopt.ipopt(model; kwargs...) print(" Second solve: "); @time sol = NLPModelsIpopt.ipopt(model; kwargs...) - println(" sol.status = $(sol.status)\n") + println(" sol.status = $(sol.status) objective (NLP) = $(sol.objective) \n") + # Test that the solver found an optimal solution - if f == :truck_trailer || - f == :space_shuttle - @test (sol.status == :infeasible) || (sol.status == :max_iter) - else - @test sol.status == :first_order + success = (sol.status == :first_order || sol.status == :acceptable) + if success + @test success + print("OptimalControl : $f converged : \033[1;32mTest Passed\033[0m\n") + else + @test success broken=true + print("OptimalControl : $f converged : \033[1;33mTest Broken\033[0m\n") + global list_of_problems_final + list_of_problems_final = setdiff(list_of_problems_final, [f]) end end end + + println() + println("\033[1m#########################################\033[0m") + println("\033[1m### END TEST CONVERGED OptimalControl ###\033[0m") + println("\033[1m#########################################\033[0m") + println() + end diff --git a/test/test_aqua.jl b/test/test_aqua.jl index 4f004e27..e430f32a 100644 --- a/test/test_aqua.jl +++ b/test/test_aqua.jl @@ -1,4 +1,10 @@ function test_aqua() + + println() + println("\033[1m###########################\033[0m") + println("\033[1m######## TEST AQUA ########\033[0m") + println("\033[1m###########################\033[0m") + println() @testset "Aqua.jl" begin Aqua.test_all( OptimalControlProblems; @@ -10,4 +16,10 @@ function test_aqua() # do not warn about ambiguities in dependencies Aqua.test_ambiguities(OptimalControlProblems) end + + println() + println("\033[1m###########################\033[0m") + println("\033[1m##### END TEST AQUA #######\033[0m") + println("\033[1m###########################\033[0m") + println() end diff --git a/test/test_custom.jl b/test/test_custom.jl new file mode 100644 index 00000000..3474e046 --- /dev/null +++ b/test/test_custom.jl @@ -0,0 +1,25 @@ +#using CTBase +using Ipopt +using JuMP +using NLPModelsIpopt +using OptimalControl +using OptimalControlProblems +using Test +include("utils.jl") + +# Parameters for the solvers +const tol = 1e-8 +const mu_strategy = "adaptive" +const sb = "yes" +const constr_viol_tol = 1e-8 +const max_iter = 1000 +const max_wall_time = 500.0 + +list_of_problems=[:electric_vehicle] +verbose = true + +include("test_JuMP.jl") +test_JuMP() + +include("test_OptimalControl.jl") +test_OptimalControl() \ No newline at end of file diff --git a/test/test_quick.jl b/test/test_quick.jl new file mode 100644 index 00000000..46be96cd --- /dev/null +++ b/test/test_quick.jl @@ -0,0 +1,50 @@ +using Printf + +function test_quick() + + # Comparison Parameters + ε = 1e-2 + + # options for solvers + kwargs = Dict( + :print_level => 0, + :tol => tol, + :mu_strategy => mu_strategy, + :sb => sb, + :max_iter => max_iter, + :max_wall_time => max_wall_time, + ) + + for f in list_of_problems + nh = OptimalControlProblems.metadata[f][:nh] + print("$f ") + ########## OptimalControl ########## + docp, OC_model = OptimalControlProblems.eval(f)(OptimalControlBackend()) + nlp_sol = NLPModelsIpopt.ipopt(OC_model; kwargs...) + sol = build_OCP_solution(docp; primal=nlp_sol.solution, dual=nlp_sol.multipliers) + obj_oc = objective(sol) + + ############### JuMP ############### + JuMP_model = OptimalControlProblems.eval(f)(JuMPBackend()) + set_optimizer(JuMP_model, Ipopt.Optimizer) + set_silent(JuMP_model) + set_optimizer_attribute(JuMP_model, "tol", tol) + set_optimizer_attribute(JuMP_model, "max_iter", max_iter) + set_optimizer_attribute(JuMP_model, "mu_strategy", mu_strategy) + set_optimizer_attribute(JuMP_model, "linear_solver", "mumps") + set_optimizer_attribute(JuMP_model, "max_wall_time", max_wall_time) + set_optimizer_attribute(JuMP_model, "sb", sb) + optimize!(JuMP_model) + obj_jmp = objective_value(JuMP_model) + + # objective relative error + dist_obj = abs(obj_oc - obj_jmp) / (obj_oc + obj_jmp) / 2 + if dist_obj < ε + @printf("Objective rel error %5.2g \033[1;32mTest Passed\033[0m\n", dist_obj) + else + @printf("Objective rel error %5.2g \033[1;33mTest Broken\033[0m JuMP: %5.2g vs OC: %5.2g\n", dist_obj, obj_jmp, obj_oc) + end + + end + +end \ No newline at end of file diff --git a/test/utils.jl b/test/utils.jl new file mode 100644 index 00000000..1d6fa0a0 --- /dev/null +++ b/test/utils.jl @@ -0,0 +1,47 @@ +function costateInterpolation(p, t) + nx = size(p[1])[1] + n_h = length(t)[1] + res = [zeros(Float64, nx) for _ in 1:n_h] + for j in 1:nx + pj = Interpolations.linear_interpolation(t[1:end-1], [p[i][j] for i in 1:n_h-1], extrapolation_bc=Interpolations.Line()) + f = t -> pj(t) + for i in 1:n_h + res[i][j] = f(t[i]) + end + end + return res +end + +function prettytime(t) + if t < 1e3 + value, units = t, "ns" + elseif t < 1e6 + value, units = t / 1e3, "μs" + elseif t < 1e9 + value, units = t / 1e6, "ms" + else + value, units = t / 1e9, "s" + end + return string(value, " ", units) +end + +### compute norm_Lp for u : R → R + +function norm_Lp(u, p, dt) + if p == Inf + if isa(u[1], Number) + return maximum(abs, u) + else + return maximum([maximum(abs, ui) for ui in u]) + end + end + nu = length(u) + if nu < 2 + return (sum(abs.(u[1]) .^ p) * dt)^(1 / p) + end + s = 0.0 + for i in 1:nu-1 + s += 0.5 * (sum(abs.(u[i]) .^ p) + sum(abs.(u[i+1]) .^ p)) + end + return (s * dt)^(1 / p) +end \ No newline at end of file