Update Optimal Control Example Documentation (#380)

* Update Fishing example docs

* Adjust unit test of fishing example to test output values within tolerance

* Update Hanging Chain example docs

* Update Jennings example docs

* Added unit tests to examples (excluding pandemic control)

* Update tolerance to 1E-06

* Minor update

---------

Co-authored-by: pulsipher <[email protected]>

Author: wenwen0231

docs/src/examples/Optimal Control/Fishing.jl

@@ -1,8 +1,9 @@
 # # Fishing Optimal Control
 # We will solve an optimal control problem to maximize profit
-# constrained by fish population
+# constrained by fish population.
+
 # # Problem Statement and Model
-# In this system out input is the rate of fishing ``u``, and the profit
+# In this system, our input is the fishing rate ``u``, and the profit
 # ``J`` will be defined in the objective to maximize.
 # Our profit objective is represented as:
 # ```math
@@ -16,51 +17,76 @@
 # &&&&&J(0) = 0 \\
 # \end{aligned}
 # ```
-# # Model Definition
+
+# # Modeling in InfiniteOpt
 # First we must import ``InfiniteOpt`` and other packages.
 using InfiniteOpt, Ipopt, Plots;
-# Next we specify an array of initial conditions as well as problem variables.
+
+# Next we specify our initial conditions and problem variables.
 x0 = 70
 E, c, r, k, Umax = 1, 17.5, 0.71, 80.5, 20;
-# We initialize the infinite model and opt to use the Ipopt solver
+
+# # Model Initialization
+# We initialize the infinite model with [`InfiniteModel`](@ref) and select Ipopt as our optimizer that will be used to solve it.
 m = InfiniteModel(Ipopt.Optimizer);
-# Now let's specify variables. ``u`` is as our fishing rate.
-# ``x`` will be used to model the fish population in response to ``u``
-# the infinite parameter ``t`` that will span over 10 years.
-@infinite_parameter(m, t in [0,10],num_supports=100)
+
+# # Infinite Parameter Definition
+# We now define the infinite parameter ``t \in [0, 10]`` to represent time over a 10 year period. We'll also specify 100 equidistant time points.
+@infinite_parameter(m, t in [0, 10], num_supports=100)
+
+# # Infinite Variable Definition
+# Now that we have our infinite parameter defined, let's specify our infinite variables:
+# - ``1 \leq x(t)`` : fish population at time ``t``
+# - ``0 \leq u(t) \leq 1`` : fishing rate
+# - ``J(t)`` : profit over time
 @variable(m, 1 <= x, Infinite(t))
 @variable(m, 0 <= u <= 1, Infinite(t));
-# ``J`` represents profit over time.
 @variable(m, J, Infinite(t));
-# Specifying the objective to maximize profit ``J``:
+
+# # Objective Definition
+# Now we add the objective using `@objective` to maximize profit ``J`` at the end of the 10 year period:
 @objective(m, Max, J(10));
-# Define the ODEs which serve as our system model.
+
+# # Constraint Definition
+# The last step is to add our constraints. First, define the ODEs which serve as our system model.
 @constraint(m, ∂(J,t) == (E-c/x) * u * Umax)
 @constraint(m, ∂(x,t) == r * x *(1 - x/k) - u*Umax);
-# Set our initial conditions.
+# We also set our initial conditions for ``x`` and ``J``.
 @constraint(m, x(0) == x0)
 @constraint(m, J(0) == 0);
+
 # # Problem Solution
-# Optimize the model:
+# Now we're ready to solve! We can solve the model by invoking `optimize!`:
 optimize!(m)
-# Extract the results.
+
+# # Extract and Plot the Results
+# Now we can extract the optimal results and plot them to visualize how the fish population compares to the profit. Note that the values of infinite variables are
+# returned as arrays corresponding to how the supports were used to discretize our model.
+# We can use the `value` function to extract the values of the infinite variables.
 ts = value(t)
 u_opt = value(u)
 x_opt = value(x)
 J_opt = value(J);
 
+# Create the plot for the fish population and profit over time.
 p1 = plot(ts, [x_opt, J_opt] ,
     label=["Fish Pop" "Profit"], 
-    title="State Variables")
+    title="State Variables");
+# Create the plot for the fishing rate over time.
 p2 = plot(ts, u_opt,
-    label = "Rate",
-    title = "Rate vs Time")
-plot(p1,p2 ,layout=(2,1), size=(800,600));
+    label = "Fishing Rate",
+    title = "Fishing Rate vs Time");
+# Visualize the two plots on one figure.
+plot(p1,p2 ,layout=(2,1), size=(800,600))
 
 # ### Maintenance Tests
 # These are here to ensure this example stays up to date. 
 using Test
+tol = 1E-6
 @test termination_status(m) == MOI.LOCALLY_SOLVED
 @test has_values(m)
 @test u_opt isa Vector{<:Real}
-@test J_opt isa Vector{<:Real}
+@test J_opt isa Vector{<:Real}
+@test isapprox(u_opt[end], 1.0000000088945395, atol=tol)
+@test isapprox(x_opt[end], 31.441105707837544, atol=tol)
+@test isapprox(J_opt[end], 106.80870543718251, atol=tol)

docs/src/examples/Optimal Control/Hanging Chain.jl

@@ -25,53 +25,78 @@
 # \end{aligned}
 # ```
 
-# # Model Definition
+# # Modeling in InfiniteOpt 
 # First we must import ``InfiniteOpt`` and other packages.
 using InfiniteOpt, Ipopt, Plots;
-# Next we specify an array of initial conditions as well as problem variables.
+
+# Next we specify our initial conditions and problem variables.
 a, b, L = 1, 3, 4
 x0 = [a, 0, 0]
 xtf = [b, NaN, L];
-# We initialize the infinite model and opt to use the Ipopt solver
+
+# # Model Definition
+# We initialize the infinite model with [`InfiniteModel`](@ref) and select Ipopt as our optimizer that will be used to solve it.
 m = InfiniteModel(Ipopt.Optimizer);
-# t is specified as ``\ t \in [0,1]``. The bounds are arbitrary for this problem:
+
+# # Infinite Parameter Definition
+# We define the infinite parameter t as ``\ t \in [0,1]``. The bounds are arbitrary for this problem:
 @infinite_parameter(m, t in [0,1], num_supports = 100);
-# Now let's specify variables. ``u`` is our controller variable.
+
+# # Infinite Variable Definition
+# Now let's specify variables. We define ``x`` as a vector containing the state variables for position, energy and chain length. ``u`` is our control variable.
 @variable(m, 0 <= x[1:3] <= 10, Infinite(t))
 @variable(m, -10 <= u <= 20, Infinite(t));
-# Specifying the objective to minimize kinetic energy at the final time:
+
+# # Objective Definition
+# We specify the objective using `@objective` to minimize potential energy at the final time:
 @objective(m, Min, x[2](1));
-# Define the ODEs which serve as our system model.
+
+# # Constraint Definition
+# The last step is to add our constraints. First, define the ODEs which serve as our system model.
 @constraint(m, ∂(x[1],t) == u)
 @constraint(m, ∂(x[2],t) == x[1] * (1 + u^2)^(1/2))
 @constraint(m, ∂(x[3],t) == (1 + u^2)^(1/2));
-# Set our inital and final conditions.
+
+# We also set our initial and final conditions for ``x``.
 @constraint(m, [i in 1:3], x[i](0) == x0[i])
 @constraint(m, x[1](1) == xtf[1])
 @constraint(m, x[3](1) == xtf[3]);
+
 # # Problem Solution
-# Optimize the model:
+# Now we can solve the model with `optimize!`:
 optimize!(m)
-# Extract the results.
+
+# # Extract and Plot the Results
+# Extract the results using `value`. Note that they are returned as arrays corresponding to the supports used to discretize the model.
 ts = value(t)
 u_opt = value(u)
 x1_opt = value(x[1])
 x2_opt = value(x[2])
-x3_opt = value(x[3])
+x3_opt = value(x[3]);
+
+# We can also check the objective value:
 @show(objective_value(m))
+
+# Create the plot for the state variables and input over time.
 p1 = plot(ts, [x1_opt, x2_opt, x3_opt], 
+    
     label=["x1" "x2" "x3"], 
     title="State Variables")
 
 p2 = plot(ts, u_opt, 
     label="u(t)", 
-    title="Input")
-plot(p1, p2, layout=(2,1), size=(800,600));
+    title="Input");
+
+# Visualize the two plots on one figure.
+plot(p1, p2, layout=(2,1), size=(800,600))
 
 # ### Maintenance Tests
 # These are here to ensure this example stays up to date. 
 using Test
+tol = 1E-6
 @test termination_status(m) == MOI.LOCALLY_SOLVED
 @test has_values(m)
 @test u_opt isa Vector{<:Real}
-@test x1_opt isa Vector{<:Real}
+@test x1_opt isa Vector{<:Real}
+@test isapprox(objective_value(m), 5.127030122851338, atol=tol)
+@test isapprox(u_opt[end], 7.20355021172144, atol=tol)

docs/src/examples/Optimal Control/Jennings.jl

@@ -27,51 +27,69 @@
 # \end{aligned}
 # ```
 
-# # Model Definition
-
+# # Modeling in InfiniteOpt
 # First we must import ``InfiniteOpt`` and other packages.
 using InfiniteOpt, Ipopt, Plots;
+
 # Next we specify an array of initial conditions.
-x0 = [π/2, 4, 0]; # x(0) for x1,x2,x3
-# We initialize the infinite model and opt to use the Ipopt solver
+x0 = [π/2, 4, 0]; # x(0) for x1, x2, x3
+
+# # Model Definition
+# We initialize the infinite model with [`InfiniteModel`](@ref), opting to use the Ipopt solver.
 m = InfiniteModel(Ipopt.Optimizer);
+
+# # Infinite Parameter Definition
 # Recall t is specified as ``\ t \in [0,1]``:
-@infinite_parameter(m, t in [0,1],num_supports= 100)
-# Now let's specify descision variables. Notice that ``t_f`` is
+@infinite_parameter(m, t in [0,1],num_supports= 100);
+
+# # Variable Definition
+# Now let's specify decision variables. Notice that ``t_f`` is
 # not a function of time and is a singular value.
 @variable(m, x[1:3], Infinite(t))
 @variable(m, -2 <= u <= 2, Infinite(t))
 @variable(m, 0.1 <= tf);
-# Specifying the objective to minimize final time:
+
+# # Objective Definition
+# Now we add the objective using `@objective` to minimize final time:
 @objective(m, Min, tf);
-# Define the ODEs which serve as our system model.
+
+# # Constraint Definition
+# The last step is to add our constraints. First, define the ODEs which serve as our system model.
 @constraint(m, ∂(x[1],t) == tf*u)
 @constraint(m, ∂(x[2],t) == tf*cos(x[1]))
 @constraint(m, ∂(x[3],t) == tf*sin(x[1]));
-# Set our inital and final conditions.
+
+# Set our inital and final conditions for ``x``.
 @constraint(m, [i in 1:3], x[i](0) == x0[i])
 @constraint(m, x[2](1) <=0)
 @constraint(m, x[3](1) <= 1e-1);
+
 # # Problem Solution
-# Optimize the model:
+# Now everything is ready for solving! We can solve the model with `@optimize!`:
 optimize!(m)
-# Extract the results. Notice that we multiply by ``t_f``
+
+# # Extract and Plot the Results
+# We can extract the results as arrays using the `value` function. Notice that we multiply by ``t_f``
 # to scale our time.
 ts = value(t)*value(tf)
 u_opt = value(u)
 x1_opt = value(x[1])
 x2_opt = value(x[2])
 x3_opt = value(x[3]);
+
 # Plot the results
 plot(ts, u_opt, label = "u(t)", linecolor = :black, linestyle = :dash)
 plot!(ts, x1_opt, linecolor = :blue, linealpha = 0.4, label = "x1")
 plot!(ts, x2_opt, linecolor = :green, linealpha = 0.4, label = "x2")
-plot!(ts, x3_opt, linecolor = :red, linealpha = 0.4, label = "x3");
+plot!(ts, x3_opt, linecolor = :red, linealpha = 0.4, label = "x3")
 
 # ### Maintenance Tests
 # These are here to ensure this example stays up to date. 
 using Test
+tol = 1E-6
 @test termination_status(m) == MOI.LOCALLY_SOLVED
 @test has_values(m)
 @test u_opt isa Vector{<:Real}
-@test x1_opt isa Vector{<:Real}
+@test x1_opt isa Vector{<:Real}
+@test isapprox(x1_opt[end], 3.2501431326448293, atol=tol)
+@test isapprox(objective_value(m), 4.284564834847627, atol=tol)

docs/src/examples/Optimal Control/consumption_savings.jl

@@ -108,7 +108,10 @@ plot!(ts[ix], BB, color = 4, linestyle=:dash, lab = "B: wealth balance, closed f
 # ### Maintenance Tests
 # These are here to ensure this example stays up to date. 
 using Test
+tol = 1E-6
 @test termination_status(m) == MOI.LOCALLY_SOLVED
 @test has_values(m)
 @test B_opt isa Vector{<:Real}
 @test c_opt isa Vector{<:Real}
+@test isapprox(opt_obj, -67025.62174598589, atol=tol)
+@test isapprox(c_opt[end], -52.00096266262752, atol=tol)

docs/src/examples/Optimal Control/hovercraft.jl

@@ -76,6 +76,11 @@ ylabel!("x_2")
 # ### Maintenance Tests
 # These are here to ensure this example stays up to date. 
 using Test
+tol = 1E-6
 @test termination_status(m) == MOI.LOCALLY_SOLVED
 @test has_values(m)
 @test x_opt isa Vector{<:Vector{<:Real}}
+@test isapprox(objective_value(m), 0.043685293177035435, atol=tol)
+@test isapprox(value(u[1])[end], -0.010503853944039986, atol=tol)
+@test isapprox(value(u[2])[end], 0.005456780217220367, atol=tol)
+