About the "horizon" parameter in the solve_ocp() method #781

zhengzl18 · 2022-10-17T13:04:26Z

zhengzl18
Oct 17, 2022

Hi,
I'm running the example code of the optimal module: https://python-control.readthedocs.io/en/0.9.2/optimal.html#example
When calling the solve_ocp method, a parameter called "horizon" is required:

horizon = np.linspace(0, Tf, 3, endpoint=True)
result = opt.solve_ocp(
    vehicle, horizon, x0, traj_cost, constraints,
    terminal_cost=term_cost, initial_guess=u0)

In this setting, it works just fine.

But when I tried to change the horizon into something like (changing nothing but the number of points in the linspace):

horizon = np.linspace(0, Tf, 4, endpoint=True)

the result seemed to be extremely strange, with the inputs u always remaining the same:

Summary statistics:
* Cost function calls: 9
* Constraint calls: 20
* System simulations: 27
* Final cost: 16000.0
    cost: 16000.0
     fun: 16000.0
  inputs: array([[10., 10., 10., 10.],
       [ 0.,  0.,  0.,  0.]])
     jac: array([ 0.00000000e+00,  2.44140625e-04,  1.22070312e-04,  0.00000000e+00,
       -3.62099047e+06, -6.84054738e+06, -2.12256850e+06, -7.49211857e+05])
 message: 'Optimization terminated successfully'
    nfev: 9
     nit: 5
    njev: 1
 problem: <control.optimal.OptimalControlProblem object at 0x7ff0cb281df0>
  states: array([[  0.        ,  33.33333333,  66.66666667, 100.        ],
       [ -2.        ,  -2.        ,  -2.        ,  -2.        ],
       [  0.        ,   0.        ,   0.        ,   0.        ]])
  status: 0
 success: True
    time: array([ 0.        ,  3.33333333,  6.66666667, 10.        ])
       x: array([10., 10., 10., 10.,  0.,  0.,  0.,  0.])

It's the case for all num except 3 in the np.linspace(0, Tf, num).

I'm confused by this results and wonder what's the actual use of 'horizon'? (I do notice from the comments that it's a "List of times at which the optimal input should be computed", but I still can't understand)

Hope for any help.

Thanks!

Answered by murrayrm

Nov 26, 2022

This problem is addressed in PR #799 by switching to using collocation as the default method for solving optimal control problems. This is much more numerically stable and provides good solutions to the example problem from the documentation (and without the issues are the solution only working for a specific number of time points).

View full answer

petterogren · 2022-10-27T11:56:56Z

petterogren
Oct 27, 2022

I am struggling with the same problem.
I think the 3 is the discretization mesh in time, thus 3 is indeed a very coarse discretization.
Furthermore, checking the variable "result" returned by the solver the "result.message" is actually "Optimization terminated successfully".
My best guess is that the tolerance in the optimization solver is too large, but I am not able to change it. I tried adding tol=0.01 but this parameter is not accepted.
Please fix this!

0 replies

murrayrm · 2022-10-27T14:33:28Z

murrayrm
Oct 27, 2022
Maintainer

I've got this on my list to take a look at, but haven't had a chance. One quick thing is that if you want to change the parameters sent to the optimizer, use the minimize_kwargs argument to solve_ocp, which lets you pass arguments to the scipy.optimize.minimize function.

1 reply

murrayrm Nov 20, 2022
Maintainer

A bit more information: the problem seems to be with the default method used by scipy.optimize.minimize, which is 'SLSQP'. Near as can tell, given the initial guess (which is to drive in a straight line), it can't figure out what to do when there are more than 3 points in the horizon. You can see that something is wrong because the value of the cost function is very high (result.cost = 16000.0 in the output above).

If you switch the solver to trust-constr, then things work OK:

horizon = np.linspace(0, Tf, 5, endpoint=True)
result = opt.solve_ocp(
    vehicle, horizon, x0, traj_cost, constraints,
    terminal_cost=term_cost, initial_guess=u0, minimize_method='trust-constr')

The (abridged) result is now

    execution_time: 226.1267969608307
               fun: 0.007588438759261349
              grad: array([-7.73912826e-03, -2.26881260e-03, -2.45970453e-03, -1.93155936e-02,
        1.53641058e-01,  2.76074462e+00,  4.69169176e+00,  3.75040997e+00,
       -2.55960229e+04,  1.59332017e+00])
            inputs: array([[ 1.00001233e+01,  1.00030576e+01,  1.00030819e+01,
         1.00008577e+01,  1.00651172e+01],
       [ 6.58643630e-03,  2.78602633e-03,  2.28824265e-03,
        -4.30400352e-03, -7.74700114e-03]])
               jac: [array([[ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -0., -0.],
       [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0., -0., -0.],
       [ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0., -0., -0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  0., -0., -0.],
       [ 0.,  0.,  1.,  0.,  0.,  0.,  0.,  0., -0., -0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1., -0., -0.],
       [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0., -0., -0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1., -0.],
       [ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0., -0., -0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -0.,  1.]])]
           problem: <control.optimal.OptimalControlProblem object at 0x10fbb6980>
            states: array([[ 0.00000000e+00,  2.49962868e+01,  4.99732645e+01,
         7.49317792e+01,  9.99999628e+01],
       [-2.00000000e+00, -1.44955389e+00, -2.03455865e-01,
         1.29540778e+00,  1.99998404e+00],
       [ 0.00000000e+00,  3.87210391e-02,  5.99980593e-02,
         5.10708705e-02,  3.22629524e-04]])
            status: 2
           success: True
              time: array([ 0. ,  2.5,  5. ,  7.5, 10. ])

The cost is now result.cost = 0.007588438759261349. Note, however, that it takes a lot longer to get the result (almost four minutes) and also you get a warning about the Hessian being approximately zero (which may be why the other solver was not working well).

I'll keep looking at this, but the quick (but slow -:) fix is to use a different optimization method.

BTW, the much faster way to solve this problem is to exploit the fact that the system is differentially flat: https://python-control.readthedocs.io/en/0.9.2/kincar-flatsys.html. Approach #3 in the example script is solving the same problem and returns a result instantaneously.

murrayrm · 2022-11-26T00:59:52Z

murrayrm
Nov 26, 2022
Maintainer

This problem is addressed in PR #799 by switching to using collocation as the default method for solving optimal control problems. This is much more numerically stable and provides good solutions to the example problem from the documentation (and without the issues are the solution only working for a specific number of time points).

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

About the "horizon" parameter in the solve_ocp() method #781

{{title}}

Replies: 3 comments 1 reply

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

Select a reply

About the "horizon" parameter in the solve_ocp() method #781

zhengzl18 Oct 17, 2022

Replies: 3 comments · 1 reply

petterogren Oct 27, 2022

murrayrm Oct 27, 2022 Maintainer

murrayrm Nov 20, 2022 Maintainer

murrayrm Nov 26, 2022 Maintainer

zhengzl18
Oct 17, 2022

Replies: 3 comments 1 reply

petterogren
Oct 27, 2022

murrayrm
Oct 27, 2022
Maintainer

murrayrm Nov 20, 2022
Maintainer

murrayrm
Nov 26, 2022
Maintainer