Issue following find_eqpt documentation #1053

corbinklett · 2024-10-30T21:03:03Z

corbinklett
Oct 30, 2024

I'm trying to use find_eqpt in order to trim an aircraft. I need to solve for the trim throttle input and the angle of attack.

There are a couple of ways that I can't figure out how to use find_eqpt, commented below. How do I specify that d/dt of only some of the states (u, w, theta, and q) matters in this calculation, given that other states are specified as fixed?

Also, it doesn't seem like find_eqpt handles bounds or constraints, such that you can constrain a quaternion to norm 1 (this example just uses Euler angles).

# Fix output at:
Vinf = 30
beta = 0
gamma = 0

# Find this output:
alpha = 0 # guess

# Find this input:
throttle = 0.5 # guess

# State guess
x = 0
y = 0
z = 0
u = Vinf
v = 0
w = 0
phi = 0    # Roll angle
theta = 0
psi = 0    # Yaw angle
p = 0
q = 0
r = 0

x0 = np.array([x, y, z, u, v, w, phi, theta, psi, p, q, r])
u0 = np.array([throttle])

trim_state, trim_input, yeq, result = ct.find_eqpt(
  aircraft_6dof_model, 
   x0, 
   u0=u0, 
   y0=[Vinf, alpha, beta, gamma],
   iy=[0, 2, 3],           			# Constrain Vinf, beta, gamma
   ix=[4, 6, 8, 9, 10, 11],   			# Constrain v, phi, psi, p, q, r
   # dx0=[0,0,0,0,0,0,0,0,0,0,0,0],	# How to specify that xdot, ydot, and zdot don't matter?
   idx=[3, 5, 7, 10], 			        # Only care about d/dt of u, w, theta, and q	
   return_y = True,
   return_result = True
)

Depending on how I formulate the arguments, I get an input/output mismatch error, or the find_eqpt function runs infinitely.

Answered by corbinklett

Nov 3, 2024

The issue is that the problem I desire to solve is over-constrained. I get the intended result if I change the line in find_eqpt that calls root to result = root(rootfun, z0, method='lm'), adding the "least squares method" designator to the args.

Does it make sense for find_eqpt to do this if num_constraints > num_freedoms?

View full answer

murrayrm · 2024-11-02T18:32:30Z

murrayrm
Nov 2, 2024
Maintainer

What you have looks correct. If you look inside the code at the definition of rootfun, you will see that the call to the root finding code has the following (extracted) lines:

dx = sys._rhs(t, x, u) - dx0
dy = sys._out(t, x, u) - y0
return np.concatenate((dx[deriv_vars], dy[output_vars]), axis=0)

The deriv_vars varaible is ix and the output_vars variable is iy, so only those variables will be constrained to be zero.

I did a run of your above using

aircraft_6dof_model = ct.rss(states=12, inputs=1, outputs=4)

and it ran OK.

If you post what you are using for aircraft_6dof_model, I can see if that is causing the problem.

In terms of including bounds or constraints: that is not supported, but I'm not sure it would help. The find_eqpt will only work if there is a unique (local) equilibrium point and so if you can't find one without the constraints, you won't find one by imposing the constraints.

5 replies

murrayrm Nov 3, 2024
Maintainer

I put your code into the latest main branch and it returned almost immediately with the following result:

In [2]: yeq
Out[2]: array([30.00017353,  0.14058246,  0.        , -0.04797305])

This doesn't quite match your hand calculation (which has 0 for the final entry), but it didn't hang and the first three entries match the final calculation that you show:

Out[1]: array([ 3.00000000e+01,  1.40499005e-01,  0.00000000e+00, -2.05860969e-17])

Can you post the version of Python, python-control, numpy, and scipy that you are using. For me, they are:

Python: 3.12.2
python-control: 0.10.1.dev32+gdbc998de
numpy: 1.26.4
scipy: 1.13.0

Also, when I look at trim_state, I get:

In [5]: trim_state
Out[5]: 
array([-6.93849578e+18, -2.68847397e+17,  0.00000000e+00,  2.97042083e+01,
        0.00000000e+00,  4.20361982e+00,  0.00000000e+00,  9.26094047e-02,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00])

The large numbers in x and y may be what is causing the problem. As you note above, those variables are "free" and so the root finding algorithm in SciPy may be having group converging.

corbinklett Nov 3, 2024
Author

The x, y large numbers may be why the result variable returns:

>>> result
 message: The iteration is not making good progress, as measured by the 
            improvement from the last five Jacobian evaluations.
 success: False
  status: 4
     fun: [ 2.526e-04 -2.221e-04  0.000e+00  0.000e+00  1.672e-04
            0.000e+00 -4.750e-02]
       x: [ 6.410e+17  7.337e+17  0.000e+00  2.970e+01  4.203e+00
            9.308e-02  2.180e-01]
  method: hybr
    nfev: 123
    fjac: [[ 1.000e+00  0.000e+00 ...  0.000e+00  0.000e+00]
           [ 0.000e+00  1.000e+00 ...  0.000e+00  0.000e+00]
           ...
           [ 0.000e+00  0.000e+00 ...  0.000e+00 -5.049e-17]
           [ 0.000e+00  0.000e+00 ...  1.000e+00  0.000e+00]]
       r: [-0.000e+00  0.000e+00 ...  0.000e+00 -0.000e+00]
     qtf: [ 2.526e-04 -2.221e-04  0.000e+00  5.682e-05 -4.750e-02
           -1.417e-18  0.000e+00]

I get the same exact results that you do, but it's still an incorrect trim result (returns a flight path angle of -2.72 degrees, but the find_eqpt y0 and iy arguments should constrain it to zero).

corbinklett Nov 3, 2024
Author

Will dig into why I can't set ix=[0, 1, 2, 4, 6, 8, 9, 10, 11] in order to fix x, y, and z, at 0.

corbinklett Nov 3, 2024
Author

The issue is that the problem I desire to solve is over-constrained. I get the intended result if I change the line in find_eqpt that calls root to result = root(rootfun, z0, method='lm'), adding the "least squares method" designator to the args.

Does it make sense for find_eqpt to do this if num_constraints > num_freedoms?

Answer selected by corbinklett

murrayrm · 2024-11-03T23:48:59Z

murrayrm
Nov 3, 2024
Maintainer

One possibility, which we do in other places, is to create a parameter root_method that let's you change the method passed to the root function. We could also include root_kwargs as a way to pass additional arguments. See ct.OptimalControlProblem for an example of this.

If we add root_method and set the default to None, then we could also add some logic to check the number of constraints versus degrees of freedom, with the default solver set to 'lm' if num_constraints > num_freedoms. Or just add something to the documentation indicating you should change the method in that case.

It still seems like setting ix appropriately should fix the problem. I'll try looking into that further when I get a chance.

1 reply

murrayrm Nov 9, 2024
Maintainer

PR #1054 adds a root_method keyword parameter that can be used to fix the issue identified above.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Issue following find_eqpt documentation #1053

{{title}}

Replies: 2 comments 6 replies

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

Issue following find_eqpt documentation #1053

corbinklett Oct 30, 2024

Replies: 2 comments · 6 replies

murrayrm Nov 2, 2024 Maintainer

murrayrm Nov 3, 2024 Maintainer

corbinklett Nov 3, 2024 Author

corbinklett Nov 3, 2024 Author

corbinklett Nov 3, 2024 Author

murrayrm Nov 3, 2024 Maintainer

murrayrm Nov 9, 2024 Maintainer

corbinklett
Oct 30, 2024

Replies: 2 comments 6 replies

murrayrm
Nov 2, 2024
Maintainer

murrayrm Nov 3, 2024
Maintainer

corbinklett Nov 3, 2024
Author

corbinklett Nov 3, 2024
Author

corbinklett Nov 3, 2024
Author

murrayrm
Nov 3, 2024
Maintainer

murrayrm Nov 9, 2024
Maintainer