Tutorial: Bilinear Control with Multiple Drives
This tutorial demonstrates a more complex bilinear control problem with:
- Multiple control inputs
- Control bounds
- A 3D state space
Problem Description
Control a 3D system with 2 independent control inputs.
Dynamics:
\[\dot{x} = (G_0 + u_1 G_1 + u_2 G_2) x\]
Goal: Navigate from [1, 0, 0] to [0, 0, 1] with bounded controls
using DirectTrajOpt
using NamedTrajectories
using LinearAlgebra
using Statistics
using PrintfStep 1: Define the System
3D System with Two Drives
Drift term - natural evolution
G_drift = [
0.0 1.0 0.0;
-1.0 0.0 0.0;
0.0 0.0 -0.1
]3×3 Matrix{Float64}:
0.0 1.0 0.0
-1.0 0.0 0.0
0.0 0.0 -0.1Drive terms - control influences
G_drive_1 = [
1.0 0.0 0.0;
0.0 0.0 0.0;
0.0 0.0 0.0
] # Controls first state
G_drive_2 = [
0.0 0.0 0.0;
0.0 0.0 1.0;
0.0 1.0 0.0
] # Controls coupling between states 2 and 3
G_drives = [G_drive_1, G_drive_2]2-element Vector{Matrix{Float64}}:
[1.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0]
[0.0 0.0 0.0; 0.0 0.0 1.0; 0.0 1.0 0.0]Generator function
G = u -> G_drift + sum(u .* G_drives)
println("System dynamics:")
println(" State dimension: 3")
println(" Number of controls: 2")
println(" G_drift creates oscillation in x₁-x₂ plane with decay in x₃")
println(" u₁ drives x₁ component")
println(" u₂ couples x₂ and x₃")System dynamics:
State dimension: 3
Number of controls: 2
G_drift creates oscillation in x₁-x₂ plane with decay in x₃
u₁ drives x₁ component
u₂ couples x₂ and x₃Step 2: Set Up the Problem
Time Discretization
N = 60
Δt = 0.15
total_time = N * Δt
println("\nTime discretization:")
println(" Time steps: $N")
println(" Step size: $Δt")
println(" Total time: $total_time seconds")
Time discretization:
Time steps: 60
Step size: 0.15
Total time: 9.0 secondsBoundary Conditions
x_init = [1.0, 0.0, 0.0]
x_goal = [0.0, 0.0, 1.0]
println("\nBoundary conditions:")
println(" Initial state: $x_init")
println(" Goal state: $x_goal")
Boundary conditions:
Initial state: [1.0, 0.0, 0.0]
Goal state: [0.0, 0.0, 1.0]Initial Guess
Linear interpolation for states
x_guess = hcat([x_init + (x_goal - x_init) * (t/(N-1)) for t = 0:(N-1)]...)3×60 Matrix{Float64}:
1.0 0.983051 0.966102 0.949153 … 0.0338983 0.0169492 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0169492 0.0338983 0.0508475 0.966102 0.983051 1.0Small random controls
u_guess = 0.1 * randn(2, N)2×60 Matrix{Float64}:
0.0708315 0.00386989 -0.0156475 … 0.108622 0.145722 -0.0277761
0.0594984 -0.00616439 -0.0933137 -0.127123 -0.170709 0.0240115Step 3: Create Trajectory with Bounds
Control bounds: -1 ≤ u ≤ 1 for both controls
traj = NamedTrajectory(
(x = x_guess, u = u_guess, Δt = fill(Δt, N));
timestep = :Δt,
controls = :u,
initial = (x = x_init,),
final = (x = x_goal,),
bounds = (u = 1.0,), # -1 ≤ u ≤ 1
)
println("\nTrajectory created:")
println(" Control bounds: ", traj.bounds.u)
Trajectory created:
Control bounds: ([-1.0, -1.0], [1.0, 1.0])Step 4: Define Dynamics and Objective
Dynamics integrator
integrator = BilinearIntegrator(G, :x, :u, traj)BilinearIntegrator: :x = exp(Δt G(:u)) :x (dim = 3)Objective: minimize control effort
R_u = 1.0 # control weight
obj = QuadraticRegularizer(:u, traj, R_u)
println("\nObjective: minimize ∫ ||u||² dt")
println(" Control weight: $R_u")
Objective: minimize ∫ ||u||² dt
Control weight: 1.0Step 5: Solve the Problem
prob = DirectTrajOptProblem(traj, obj, integrator)
probDirectTrajOptProblem
Trajectory
Timesteps: 60
Duration: 8.85
Knot dim: 6
Variables: x (3), u (2), Δt (1)
Controls: u, Δt
Objective: QuadraticRegularizer on :u (R = [1.0, 1.0], all)
Dynamics (1 integrators)
BilinearIntegrator: :x = exp(Δt G(:u)) :x (dim = 3)
Constraints (4 total: 2 equality, 2 bounds)
EqualityConstraint: "initial value of x"
EqualityConstraint: "final value of x"
BoundsConstraint: "bounds on u"
BoundsConstraint: "bounds on Δt"println("Solving optimization problem...")
println("="^50)
solve!(prob; max_iter = 150, verbose = false)
println("="^50)
println("Optimization complete!")
println("="^50)Solving optimization problem...
==================================================
This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.
Number of nonzeros in equality constraint Jacobian...: 2106
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 3318
Total number of variables............................: 354
variables with only lower bounds: 60
variables with lower and upper bounds: 120
variables with only upper bounds: 0
Total number of equality constraints.................: 177
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.4170545e-02 1.50e-01 1.00e+00 0.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 4.9924072e-03 1.13e-01 1.34e+00 -4.0 5.95e-01 0.0 7.04e-01 2.49e-01h 1
2 1.0853317e-02 1.06e-01 3.59e+00 -4.0 2.48e+00 0.4 1.10e-01 5.63e-02h 1
3 5.4360912e-02 9.66e-02 4.97e+01 -0.4 3.07e+00 0.9 8.86e-01 9.20e-02h 1
4 1.2004229e-01 8.86e-02 9.78e+01 -0.0 1.88e+00 1.3 1.00e+00 1.22e-01h 1
Number of Iterations....: 5
Number of objective function evaluations = 6
Number of objective gradient evaluations = 6
Number of equality constraint evaluations = 6
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 5
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 5
Total seconds in IPOPT = 6.657
EXIT: Invalid number in NLP function or derivative detected.
==================================================
Optimization complete!
==================================================Step 6: Analyze the Solution
x_sol = prob.trajectory.x
u_sol = prob.trajectory.u
times = cumsum([0.0; prob.trajectory.Δt[:]])61-element Vector{Float64}:
0.0
0.1113210124759489
0.24529345084663423
0.46840163521503053
0.5378668121780527
0.561423658104421
0.5923716246216131
1.0815780025638522
1.089230344782912
1.0954978943784481
⋮
2.1913063497035767
2.1974629098434955
2.203554025280226
2.2096202519188424
2.216268390483498
2.2230917071055085
2.2285487808683353
2.233539717058609
2.45426009171469Goal Reaching
println("\nGoal reaching:")
println(" Initial state: ", x_sol[:, 1])
println(" Final state: ", x_sol[:, end])
println(" Goal state: ", x_goal)
println(" Final error: ", norm(x_sol[:, end] - x_goal))
Goal reaching:
Initial state: [1.0, 0.0, 0.0]
Final state: [0.0, 0.0, 1.0]
Goal state: [0.0, 0.0, 1.0]
Final error: 0.0Control Statistics
println("\nControl statistics:")
for i = 1:2
u_i = u_sol[i, :]
println(" u$i:")
println(" Max magnitude: ", maximum(abs.(u_i)))
println(" Mean magnitude: ", mean(abs.(u_i)))
println(" Total norm: ", norm(u_i))
# Check bound satisfaction
if all(-1.0 .<= u_i .<= 1.0)
println(" ✓ Bounds satisfied")
else
println(" ✗ Bounds violated!")
end
end
Control statistics:
u1:
Max magnitude: 0.27353529173333735
Mean magnitude: 0.09888866566236919
Total norm: 0.9328674497430671
✓ Bounds satisfied
u2:
Max magnitude: 0.9985745674875376
Mean magnitude: 0.2052466401497687
Total norm: 2.3164120929469423
✓ Bounds satisfiedState Trajectory Analysis
println("\nState trajectory:")
println(" Max |x₁|: ", maximum(abs.(x_sol[1, :])))
println(" Max |x₂|: ", maximum(abs.(x_sol[2, :])))
println(" Max |x₃|: ", maximum(abs.(x_sol[3, :])))
State trajectory:
Max |x₁|: 1.0062079940710362
Max |x₂|: 0.6275766571499769
Max |x₃|: 1.0Step 7: Detailed Results
println("\n" * "="^50)
println("SOLUTION DETAILS")
println("="^50)
println("\nState trajectory (selected time points):")
println("Time | x₁ | x₂ | x₃")
println("-"^40)
for k in [1, 10, 20, 30, 40, 50, N]
t = times[k]
println(
@sprintf("%.2f | %7.4f | %7.4f | %7.4f", t, x_sol[1, k], x_sol[2, k], x_sol[3, k])
)
end
println("\nControl trajectory (selected time points):")
println("Time | u₁ | u₂")
println("-"^30)
for k in [1, 10, 20, 30, 40, 50, N]
t = times[k]
println(@sprintf("%.2f | %7.4f | %7.4f", t, u_sol[1, k], u_sol[2, k]))
end
==================================================
SOLUTION DETAILS
==================================================
State trajectory (selected time points):
Time | x₁ | x₂ | x₃
----------------------------------------
0.00 | 1.0000 | 0.0000 | 0.0000
1.10 | 0.7916 | -0.5176 | 0.2730
1.51 | 0.5602 | -0.3311 | 0.4563
2.05 | 0.3903 | -0.3349 | 0.6089
2.11 | 0.2816 | -0.1279 | 0.7292
2.17 | 0.1218 | -0.0716 | 0.8588
2.23 | 0.0000 | 0.0000 | 1.0000
Control trajectory (selected time points):
Time | u₁ | u₂
------------------------------
0.00 | 0.1760 | -0.0654
1.10 | -0.0730 | -0.0993
1.51 | 0.1706 | -0.2642
2.05 | -0.2207 | -0.2030
2.11 | 0.1917 | -0.0710
2.17 | -0.0162 | -0.1383
2.23 | -0.0268 | 0.0231Step 8: Verify Dynamics Satisfaction
println("\nDynamics verification:")
max_error = 0.0
for k = 1:(N-1)
global max_error
x_k = x_sol[:, k]
u_k = u_sol[:, k]
Δt_k = prob.trajectory.Δt[k]
# Predicted next state
x_k1_pred = exp(Δt_k * G(u_k)) * x_k
x_k1_actual = x_sol[:, k+1]
err = norm(x_k1_pred - x_k1_actual)
max_error = max(max_error, err)
end
println(" Maximum dynamics error: $max_error")
Dynamics verification:
Maximum dynamics error: 0.09444487985121011Comparison: Different Control Weights
println("\n" * "="^50)
println("EXPLORING DIFFERENT CONTROL WEIGHTS")
println("="^50)
==================================================
EXPLORING DIFFERENT CONTROL WEIGHTS
==================================================Try with higher control penalty
traj_high = NamedTrajectory(
(x = x_guess, u = 0.1*randn(2, N), Δt = fill(Δt, N));
timestep = :Δt,
controls = :u,
initial = (x = x_init,),
final = (x = x_goal,),
bounds = (u = 1.0,),
)
obj_high = QuadraticRegularizer(:u, traj_high, 10.0) # 10x larger weight
integrator_high = BilinearIntegrator(G, :x, :u, traj_high)
prob_high = DirectTrajOptProblem(traj_high, obj_high, integrator_high)
println("\nSolving with high control weight (R=10.0)...")
solve!(prob_high; max_iter = 150, verbose = false)
u_sol_high = prob_high.trajectory.u
println("\nComparison:")
println(" Original (R=1.0):")
println(" ||u||: ", norm(u_sol))
println(" High weight (R=10.0):")
println(" ||u||: ", norm(u_sol_high))
println(" Control reduction: ", (1 - norm(u_sol_high)/norm(u_sol)) * 100, "%")┌ Warning: Trajectory has timestep variable :Δt but no bounds on it.
│ Adding default lower bound of 0 to prevent negative timesteps.
│
│ Recommended: Add explicit bounds when creating the trajectory:
│ NamedTrajectory(...; Δt_bounds=(min, max))
│ Example:
│ NamedTrajectory(qtraj, N; Δt_bounds=(1e-3, 0.5))
│
│ Or use timesteps_all_equal=true in problem options to fix timesteps.
└ @ DirectTrajOpt.Problems ~/work/DirectTrajOpt.jl/DirectTrajOpt.jl/src/problems.jl:66
Solving with high control weight (R=10.0)...
This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.
Number of nonzeros in equality constraint Jacobian...: 2106
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 3318
Total number of variables............................: 354
variables with only lower bounds: 60
variables with lower and upper bounds: 120
variables with only upper bounds: 0
Total number of equality constraints.................: 177
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.0153459e-01 1.50e-01 1.05e+00 0.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 3.8980493e-02 1.16e-01 1.58e+00 -4.0 6.57e-01 0.0 6.77e-01 2.26e-01h 1
2 1.0627096e-01 1.07e-01 7.35e+00 -0.6 2.53e+00 0.4 1.94e-01 7.97e-02h 1
3 4.5232991e-01 1.02e-01 3.44e+01 -0.1 6.46e+00 0.9 5.77e-01 4.28e-02h 1
4 1.4784777e+00 1.08e-01 1.63e+02 0.4 2.72e+00 1.3 8.47e-01 1.23e-01h 1
5 1.5892854e+00 1.74e-01 3.81e+02 1.4 1.00e+01 0.8 6.50e-01 6.92e-02f 1
6 3.1261983e+00 1.18e-01 5.39e+02 1.1 1.41e+00 - 6.78e-01 4.20e-01h 1
7 3.5641988e+00 8.27e-02 4.90e+02 0.7 2.28e+00 - 3.18e-01 2.29e-01h 1
8 6.3399132e+00 1.23e-01 9.24e+02 1.7 8.70e+00 0.3 1.01e-01 8.45e-02h 1
9 6.5236226e+00 1.19e-01 4.50e+04 1.6 1.06e+00 3.5 1.00e+00 3.24e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.0018882e+01 3.69e-01 4.13e+04 2.9 4.05e+01 - 1.40e-01 1.63e-01f 1
11 9.5923890e+00 3.35e-01 1.40e+05 2.1 1.40e+00 3.9 7.42e-01 8.84e-02h 1
12 9.4211668e+00 3.01e-01 1.44e+05 1.4 1.50e+00 3.4 4.03e-01 1.11e-01h 1
13 9.5641474e+00 2.86e-01 1.33e+05 1.4 1.26e+01 2.9 9.84e-02 5.16e-02h 1
14 1.1823363e+01 2.32e-01 2.87e+05 2.5 1.70e+00 3.4 1.00e+00 2.16e-01h 1
15 1.2064585e+01 2.31e-01 2.76e+05 1.3 8.34e+00 2.9 1.99e-01 8.43e-03h 4
16 1.2176646e+01 2.29e-01 2.71e+05 -3.5 7.70e+00 3.3 3.73e-03 8.01e-03h 4
17 1.2370154e+01 2.27e-01 4.49e+06 2.8 4.74e+00 4.6 1.00e+00 7.44e-03h 1
18 1.2746294e+01 2.24e-01 4.15e+06 -2.8 8.42e+00 - 8.28e-02 1.28e-02h 3
19 1.3507232e+01 2.17e-01 4.12e+06 -2.8 5.31e+00 - 1.41e-03 3.05e-02h 2
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 1.4482302e+01 2.15e-01 4.63e+06 2.9 6.93e+00 4.2 9.44e-02 1.32e-02h 2
21 1.4110273e+01 2.09e-01 4.62e+06 2.9 1.61e+02 - 1.06e-03 1.14e-02h 1
22 1.3977312e+01 2.07e-01 4.72e+06 2.9 1.19e+01 - 7.20e-02 5.80e-03h 1
23 1.3034427e+01 1.97e-01 8.09e+06 2.6 4.99e+00 3.7 2.27e-01 3.94e-02h 1
24 1.4075915e+01 1.43e-01 8.15e+06 2.9 1.61e+01 5.0 5.32e-03 1.75e-02h 1
25 1.4164772e+01 1.43e-01 7.30e+06 2.9 1.23e+01 4.5 1.10e-01 1.40e-03h 1
26 1.7367561e+01 1.47e-01 7.34e+06 2.9 1.94e+01 - 4.92e-02 4.63e-02h 1
27 1.7289664e+01 1.46e-01 8.00e+06 -2.7 7.55e+01 - 1.20e-02 6.76e-03h 1
28 1.7346860e+01 1.46e-01 8.02e+06 2.9 1.43e+03 - 8.37e-04 5.76e-04h 1
29 1.7436147e+01 1.46e-01 1.04e+07 -2.7 7.88e+01 - 7.84e-03 9.79e-04h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 1.7540620e+01 1.46e-01 5.49e+247 2.9 1.44e+03 - 3.80e-04 1.01e-03h 1
31 1.7537635e+01 1.46e-01 9.62e+06 1.0 2.72e+244 - 3.63e-245 1.43e-247h 9
32 1.7625278e+01 1.46e-01 2.62e+07 -2.6 7.29e+01 - 1.10e-02 8.73e-04h 1
33 1.7751067e+01 1.45e-01 2.72e+07 -2.4 8.38e+01 - 9.12e-05 4.44e-04h 1
34 1.8035538e+01 1.45e-01 2.80e+07 2.9 5.33e+01 - 9.53e-04 1.44e-03f 1
35 3.1797139e+01 1.27e-01 6.08e+07 2.9 1.51e+01 - 5.62e-02 6.51e-02f 1
36 3.1835051e+01 1.27e-01 6.06e+07 2.9 1.88e+00 6.8 9.17e-03 3.42e-03h 1
37 3.1993981e+01 1.24e-01 5.92e+07 -2.6 2.37e+00 6.3 9.36e-02 2.12e-02h 1
38 3.2239966e+01 1.19e-01 5.64e+07 2.9 1.55e+00 6.7 8.59e-03 4.14e-02f 1
39 3.2106310e+01 1.17e-01 8.52e+116 2.9 2.13e+00 7.1 6.02e-03 1.85e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 8.3706003e+01 1.17e-01 5.54e+07 1.0 2.32e+109 7.6 4.34e-109 1.37e-106h 17
41 8.1888486e+01 1.15e-01 5.46e+07 2.9 2.81e+00 7.1 8.80e-03 1.65e-02f 1
42 8.1255123e+01 1.14e-01 5.43e+07 2.9 2.01e+00 7.5 1.58e-02 7.52e-03h 1
43 8.0127751e+01 1.13e-01 5.39e+07 2.9 5.02e+00 7.0 5.72e-03 1.02e-02h 1
44 7.5309400e+01 1.07e-01 5.06e+07 2.9 1.88e+00 7.5 3.42e-03 5.85e-02h 1
45 7.4811059e+01 1.07e-01 5.05e+07 1.8 1.36e+00 7.9 3.35e-02 8.56e-03h 1
46 7.2439707e+01 1.03e-01 4.94e+07 -2.6 2.74e+00 7.4 2.27e-03 3.00e-02f 1
47 6.9097605e+01 9.72e-02 8.52e+116 2.9 1.55e+00 7.8 3.60e-03 5.96e-02f 1
48 6.2800471e+01 9.72e-02 4.53e+07 1.0 4.64e+108 8.3 6.80e-106 1.27e-105h 11
49 5.6681809e+01 8.78e-02 4.28e+07 2.9 1.24e+00 7.8 1.06e-01 9.74e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 5.5608475e+01 8.59e-02 4.43e+07 2.9 6.78e-01 8.2 4.53e-02 2.11e-02h 1
51 5.2967776e+01 8.17e-02 4.36e+07 2.9 7.99e-01 7.7 3.82e-02 4.88e-02f 1
52 5.1276794e+01 7.87e-02 4.49e+07 2.9 5.27e-01 8.2 1.00e-01 3.71e-02h 1
53 4.8803099e+01 7.46e-02 4.33e+07 2.9 7.07e-01 7.7 1.88e-01 5.21e-02f 1
54 4.3197023e+01 5.72e-02 4.69e+07 2.9 4.73e-01 8.1 4.49e-02 2.33e-01h 1
55 4.2796462e+01 5.63e-02 4.73e+07 2.9 2.83e-01 8.5 1.36e-01 1.60e-02h 1
56 2.9595155e+01 3.34e-02 3.40e+07 2.9 3.41e-01 8.1 2.16e-01 6.15e-01f 1
57 3.1669156e+01 5.59e-02 4.67e+07 2.9 2.73e-01 7.6 1.77e-01 7.11e-01f 1
58 2.6411250e+01 6.75e-04 8.64e+06 2.9 3.60e-02 7.1 8.74e-01 1.00e+00f 1
59 2.6386888e+01 3.20e-06 1.22e+05 2.2 2.34e-03 6.6 9.89e-01 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 2.6386796e+01 5.86e-12 1.43e+03 -3.8 5.45e-06 6.2 9.90e-01 1.00e+00h 1
61 9.5912879e+00 1.18e+00 1.27e+02 1.0 4.84e+04 - 9.16e-01 8.06e-01f 1
62 9.4731182e+00 1.11e+00 1.61e+05 0.7 7.29e-01 5.7 1.52e-01 5.46e-02h 1
63 9.4694240e+00 1.11e+00 1.35e+05 0.7 6.53e-01 5.2 1.00e+00 1.67e-03h 1
64 9.2362992e+00 9.40e-01 1.11e+05 0.7 6.01e-01 4.7 2.08e-01 1.48e-01h 1
65 9.3100152e+00 8.61e-01 1.21e+05 0.7 5.23e-01 5.1 1.00e+00 8.40e-02h 1
66 9.3271288e+00 8.45e-01 2.47e+05 0.7 6.68e-01 5.6 6.14e-01 1.77e-02h 1
67 1.1064448e+01 5.00e-01 6.80e+04 0.7 6.57e-01 5.1 1.00e+00 3.97e-01h 1
68 1.1376334e+01 4.46e-01 1.37e+05 0.7 3.55e-01 5.5 1.00e+00 1.08e-01h 1
69 1.1748085e+01 3.95e-01 6.22e+04 0.7 3.38e-01 5.0 5.20e-01 1.13e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
70 1.1924991e+01 3.63e-01 1.96e+05 0.7 2.69e-01 5.5 1.00e+00 8.12e-02h 1
71 1.2473018e+01 2.89e-01 1.11e+05 0.7 2.81e-01 5.0 4.91e-01 2.02e-01h 1
72 1.2740637e+01 2.43e-01 1.09e+05 0.7 1.98e-01 5.4 1.00e+00 1.58e-01h 1
73 1.3327324e+01 1.62e-01 3.68e+04 0.7 1.84e-01 4.9 7.73e-01 3.31e-01h 1
74 1.3780033e+01 8.34e-02 2.98e+04 0.3 1.09e-01 5.4 1.00e+00 4.79e-01h 1
75 1.3919115e+01 5.92e-02 4.53e+04 -0.3 5.93e-02 4.9 1.00e+00 2.88e-01h 1
76 1.4097601e+01 2.66e-02 5.26e+04 -1.1 3.97e-02 5.3 1.00e+00 5.47e-01h 1
77 1.4110640e+01 2.39e-02 6.25e+09 -2.4 1.83e-02 4.8 1.00e+00 1.01e-01h 1
78 1.4110256e+01 2.39e-02 4.02e+04 0.3 1.98e+03 5.3 5.32e-07 2.48e-06f 1
Number of Iterations....: 79
Number of objective function evaluations = 134
Number of objective gradient evaluations = 80
Number of equality constraint evaluations = 134
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 79
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 79
Total seconds in IPOPT = 1.441
EXIT: Invalid number in NLP function or derivative detected.
Comparison:
Original (R=1.0):
||u||: 2.49719976436427
High weight (R=10.0):
||u||: 7.875755734122132
Control reduction: -215.3834885983629%Key Observations
- Multiple controls allow independent actuation of different system modes
- Bounds are strictly enforced — check that max|u| ≤ 1
- Control weight affects aggressiveness: lower weight produces larger controls that may saturate bounds, higher weight produces gentler controls
- Initial guess matters for bounded problems — random works here, but more complex problems may need better initialization
- BilinearIntegrator handles multi-input systems naturally — just provide all drive matrices in a vector
Exercises
Try these modifications:
Exercise 1: Asymmetric Bounds
Set different bounds for each control:
bounds=(u = ([-1.0, -0.5], [1.0, 2.0]),)Exercise 2: Initial Control Constraints
Start and end with zero control:
initial=(x = x_init, u = [0.0, 0.0]),
final=(x = x_goal, u = [0.0, 0.0])Exercise 3: Intermediate Waypoint
Add a waypoint constraint at t=30:
waypoint = [0.5, 0.5, 0.5]
constraint = NonlinearKnotPointConstraint(
x -> x - waypoint, :x, traj;
times=[30], equality=true
)Exercise 4: Different Goal
Try reaching [0, 1, 0] instead of [0, 0, 1]
Exercise 5: Add Minimum Time
Make Δt variable and add MinimumTimeObjective:
obj = QuadraticRegularizer(:u, traj, 1.0) +
MinimumTimeObjective(traj, 0.5)
bounds=(u = 1.0, Δt = (0.05, 0.3))Next Steps
- Minimum Time Tutorial: Optimize trajectory duration
- Smooth Controls Tutorial: Add derivative penalties for implementability
This page was generated using Literate.jl.