Shortest parametric curve in the plane

[dpo]

Hello. In my class, one of the simplest examples of a control problem that I use is to find the shortest curve $(x(t), y(t))$ in the plane between two fixed points. Of course, the solution should be a line segment. Here is my model:

t0 = 0
tf = 1
x0 = 0
xf = 1
y0 = 0
yf = 1

ocp = @def begin
    t ∈ [ t0, tf ], time
    x ∈ R², state
    u ∈ R², control

    x(t0) == [x0, y0]
    x(tf) == [xf, yf]

    ẋ(t) == [u1(t), u2(t)]

    ∫( √(u1(t)^2 + u2(t)^2) ) → min
end

I solve the model with IPOPT:

sol = solve(ocp)

This is Ipopt version 3.14.17, running with linear solver MUMPS 5.8.0.

Number of nonzeros in equality constraint Jacobian...:     3505
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      753

Total number of variables............................:     1255
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:      755
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.0000000e-01 9.00e-01 6.77e-15   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.4142136e+00 4.44e-16 1.41e-04 -11.0 1.41e+00  -4.0 1.00e+00 1.00e+00h  1
   2  1.4142136e+00 2.22e-16 1.30e-18 -11.0 2.68e-14  -4.5 1.00e+00 1.00e+00h  1

Number of Iterations....: 2

                                   (scaled)                 (unscaled)
Objective...............:   1.4142135623730947e+00    1.4142135623730947e+00
Dual infeasibility......:   1.3010426069826053e-18    1.3010426069826053e-18
Constraint violation....:   2.2204460492503131e-16    2.2204460492503131e-16
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   2.2204460492503131e-16    2.2204460492503131e-16


Number of objective function evaluations             = 3
Number of objective gradient evaluations             = 3
Number of equality constraint evaluations            = 3
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 3
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 2
Total seconds in IPOPT                               = 0.998

EXIT: Optimal Solution Found.

The final objective value makes sense. But when I plot the solution,

plot(sol)

I get

So ... not exactly a straight line. Did I do something wrong?

Thanks!

[ocots]

Actually you get a straight line:

x = state(sol)
x1 = t -> x(t)[1]
x2 = t -> x(t)[2]
t = time_grid(sol)
plot(x1.(t), x2.(t), xlabel="x₁", ylabel="x₂", legend=:none)

However, with the square root in the integrand of the objective function, I think that the problem is not well-posed. You should compute the solution of:

ocp = @def begin
    t ∈ [ t0, tf ], time
    x ∈ R², state
    u ∈ R², control

    x(t0) == [x0, y0]
    x(tf) == [xf, yf]

    ẋ(t) == [u1(t), u2(t)]

    ∫( u1(t)^2 + u2(t)^2 ) → min
end

which is equivalent to find the solution of your problem with constant control.

[ocots]

Consider the problem:

$$\min \int_{t_0}^{t_f} \| u(t) \| \mathrm{d}t, \quad \dot{x}(t) = u(t), \quad x(t_0) = A, \quad x(t_f) = B.$$

The pseudo-Hamiltonian (normal case $p^0=-1$) is

$$H(x, p, u) = p u - \| u \|$$

The costate is constant since $\partial_x H = 0$. We write $p(t) = w$. The maximisation condition gives:

$$\frac{\partial H}{\partial u}(x, p, u) = p - \frac{u}{\| u \|} = 0$$

Hence, $| p(t) | = | w | = 1$. You can check this with:

p = costate(sol)
using LinearAlgebra
norm(p(t0))

Besides, $u(t) = \lambda(t) p(t) = \lambda w$. The length of the path is (that is the cost):

$$\int_{t_0}^{t_f} \| u(t) \| \mathrm{d}t = \int_{t_0}^{t_f} \lambda(t) \mathrm{d}t.$$

Now $\lambda$ is constrained by the fact that we want to join $A$ and $B$. We have:

$$\dot{x}(t) = u(t) = \lambda(t) w$$

so

$$B - A = w \int_{t_0}^{t_f} \lambda(t) \mathrm{d}t$$

or

$$\int_{t_0}^{t_f} \lambda(t) \mathrm{d}t = \frac{\| B-A \|}{\| w \|} = \| B-A \|.$$

At the end, any $\lambda$ that satisfy this condition is a solution and the cost is $| B-A |$.

[ocots]

Another possibility is to restrict to (parameterized) paths with constant speed:

ocp = @def begin
    t ∈ [ t0, tf ], time
    q = (x₁, x₂, λ) ∈ R³, state
    u ∈ R², control

    [x₁, x₂](t0) == [x0, y0]
    [x₁, x₂](tf) == [xf, yf]

    λ(t) - (u1(t)^2 + u2(t)^2) == 0

    q̇(t) == [
        u1(t), 
        u2(t), 
        0*λ(t)      # sould be 0 but it bugs because of AD
    ] 

    ∫( √(u1(t)^2 + u2(t)^2) ) → min
end

this is equivalent to replace ∫( √(u1(t)^2 + u2(t)^2) ) by ∫( u1(t)^2 + u2(t)^2 ).

[dpo]

Thank you. Actually, a simple solution is to change the objective to

∫( 1 + u1(t)^2 + u2(t)^2 ) → min

i.e., we consider that the curve lies in a plane in $\mathbb{R}^3$. That’s obviously an equivalent problem and the solution looks great.

[ocots]

Have you tried simply:

ocp = @def begin
    t ∈ [ t0, tf ], time
    x ∈ R², state
    u ∈ R², control

    x(t0) == [x0, y0]
    x(tf) == [xf, yf]

    ẋ(t) == [u1(t), u2(t)]

    ∫( u1(t)^2 + u2(t)^2 ) → min
end

[dpo]

Not yet. I need to convince myself that the two problems are actually equivalent.

In my reply above, I copy-pasted the wrong objective. I meant

∫(√(1 + u1(t)^2 + u2(t)^2) ) → min

[ocots]

@dpo Another point of view (maybe the best) is to consider the minimum time problem, but we need to impose some constraints on the control.

t0 = 0
A = [0, 0]
B = [1, 1]

ocp = @def begin
    tf ∈ R, variable
    t ∈ [ t0, tf ], time
    x ∈ R², state
    u ∈ R², control

    tf ≥ 0.1 # to help convergence

    x(t0) == A
    x(tf) == B

    u1(t)^2 + u2(t)^2 ≤ 1

    ẋ(t) == u(t)

    tf → min
end

The final time tf solution of this problem gives the length of the path, since the solution satisfies u1(t)^2 + u2(t)^2 = 1 and so the curve is parameterized with arc length.