GeezerVison is showing github.com — Issue following find_eqpt documentation · python-control/python-control · Discussion #1053
Issue following find_eqpt documentation
#1053
Oct 30, 2024
·
2 comments
·
6 replies
Discussion options
-
Uh oh!
There was an error while loading. Please reload this page.
|
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.
|
You must be logged in to vote
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?
Replies:
2 comments
·
6 replies
Comment options
-
Uh oh!
There was an error while loading. Please reload this page.
|
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.
|
You must be logged in to vote
5 replies
Comment options
-
Uh oh!
There was an error while loading. Please reload this page.
-
Uh oh!
There was an error while loading. Please reload this page.
|
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.
|
Comment options
-
Uh oh!
There was an error while loading. Please reload this page.
|
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).
|
Comment options
-
Uh oh!
There was an error while loading. Please reload this page.
|
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.
|
Comment options
-
Uh oh!
There was an error while loading. Please reload this page.
|
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?
|
Comment options
-
Uh oh!
There was an error while loading. Please reload this page.
|
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.
|
You must be logged in to vote
1 reply
Comment options
-
Uh oh!
There was an error while loading. Please reload this page.
|
PR #1054 adds a root_method keyword parameter that can be used to fix the issue identified above.
|
|