8000 Merge pull request #1054 from murrayrm/find_operating_point-06Nov2024 · python-control/python-control@d810d5a · GitHub
[go: up one dir, main page]
Skip to content

Commit d810d5a

Browse files
murrayrmmurrayrm
authored
Merge pull request #1054 from murrayrm/find_operating_point-06Nov2024
Update find_eqpt to find_operating_point, adding root_method + better docs
2 parents cc0020a + e6fdc17 commit d810d5a

File tree

12 files changed

+323
-93
lines changed

12 files changed

+323
-93
lines changed

‎control/matlab/__init__.py

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -84,10 +84,12 @@
8484
from ..dtime import c2d
8585
from ..sisotool import sisotool
8686
from ..stochsys import lqe, dlqe
87+
from ..nlsys import find_operating_point
8788

8889
# Functions that are renamed in MATLAB
8990
pole, zero = poles, zeros
9091
freqresp = frequency_response
92+
trim = find_operating_point
9193

9294
# Import functions specific to Matlab compatibility package
9395
from .timeresp import *

‎control/nlsys.py

Lines changed: 198 additions & 71 deletions
Large diffs are not rendered by default.
Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,8 @@
3939
from .ctrlplot import ControlPlot, _add_arrows_to_line2D, _get_color, \
4040
_process_ax_keyword, _update_plot_title
4141
from .exception import ControlNotImplemented
42-
from .nlsys import NonlinearIOSystem, find_eqpt, input_output_response
42+
from .nlsys import NonlinearIOSystem, find_operating_point, \
43+
input_output_response
4344

4445
__all__ = ['phase_plane_plot', 'phase_plot', 'box_grid']
4546

@@ -853,7 +854,7 @@ def _find_equilpts(sys, points, params=None):
853854
equilpts = []
854855
for i, x0 in enumerate(points):
855856
# Look for an equilibrium point near this point
856-
xeq, ueq = find_eqpt(sys, x0, 0, params=params)
857+
xeq, ueq = find_operating_point(sys, x0, 0, params=params)
857858

858859
if xeq is None:
859860
continue # didn't find anything

‎control/tests/docstrings_test.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,7 @@
5353
control.create_estimator_iosystem: ['state_labels'], # deprecated
5454
control.bode_plot: ['sharex', 'sharey', 'margin_info'], # deprecated
5555
control.eigensys_realization: ['arg'], # quasi-positional
56+
control.find_operating_point: ['method'], # internal use
5657
}
5758

5859
# Decide on the level of verbosity (use -rP when running pytest)

‎control/tests/iosys_test.py

Lines changed: 97 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -10,10 9E88 +10,11 @@
1010

1111
import re
1212
import warnings
13-
import pytest
13+
from math import sqrt
1414

1515
import numpy as np
16-
from math import sqrt
16+
import pytest
17+
import scipy
1718

1819
import control as ct
1920

@@ -2087,6 +2088,100 @@ def test_find_eqpt(x0, ix, u0, iu, y0, iy, dx0, idx, dt, x_expect, u_expect):
20872088
np.testing.assert_allclose(np.array(ueq), u_expect, atol=1e-6)
20882089

20892090

2091+
# Test out new operating point version of find_eqpt
2092+
def test_find_operating_point():
2093+
dt = 1
2094+
sys = ct.NonlinearIOSystem(
2095+
eqpt_rhs, eqpt_out, dt=dt, states=3, inputs=2, outputs=2)
2096+
2097+
# Conditions that lead to no exact solution (from previous unit test)
2098+
x0 = 0; ix = None
2099+
u0 = [-1, 0]; iu = None
2100+
y0 = None; iy = None
2101+
dx0 = None; idx = None
2102+
2103+
# Default version: no equilibrium solution => returns None
2104+
op_point = ct.find_operating_point(
2105+
sys, x0, u0, y0, ix=ix, iu=iu, iy=iy, dx0=dx0, idx=idx)
2106+
assert op_point.states is None
2107+
assert op_point.inputs is None
2108+
assert op_point.result.success is False
2109+
2110+
# Change the method to Levenberg-Marquardt (gives nearest point)
2111+
op_point = ct.find_operating_point(
2112+
sys, x0, u0, y0, ix=ix, iu=iu, iy=iy, dx0=dx0, idx=idx,
2113+
root_method='lm')
2114+
assert op_point.states is not None
2115+
assert op_point.inputs is not None
2116+
assert op_point.result.success is True
2117+
2118+
# Make sure we get a solution if we ask for the result explicitly
2119+
op_point = ct.find_operating_point(
2120+
sys, x0, u0, y0, ix=ix, iu=iu, iy=iy, dx0=dx0, idx=idx,
2121+
return_result=True)
2122+
assert op_point.states is not None
2123+
assert op_point.inputs is not None
2124+
assert op_point.result.success is False
2125+
2126+
# Check to make sure unknown keywords are caught
2127+
with pytest.raises(TypeError, match="unrecognized keyword"):
2128+
ct.find_operating_point(sys, x0, u0, unknown=None)
2129+
2130+
2131+
def test_operating_point():
2132+
dt = 1
2133+
sys = ct.NonlinearIOSystem(
2134+
eqpt_rhs, eqpt_out, dt=dt, states=3, inputs=2, outputs=2)
2135+
2136+
# Find the operating point near the origin
2137+
op_point = ct.find_operating_point(sys, 0, 0)
2138+
2139+
# Linearize the old fashioned way
2140+
linsys_orig = ct.linearize(sys, op_point.states, op_point.inputs)
2141+
2142+
# Linearize around the operating point
2143+
linsys_oppt = ct.linearize(sys, op_point)
2144+
2145+
np.testing.assert_allclose(linsys_orig.A, linsys_oppt.A)
2146+
np.testing.assert_allclose(linsys_orig.B, linsys_oppt.B)
2147+
np.testing.assert_allclose(linsys_orig.C, linsys_oppt.C)
2148+
np.testing.assert_allclose(linsys_orig.D, linsys_oppt.D)
2149+
2150+
# Call find_operating_point with method and keyword arguments
2151+
op_point = ct.find_operating_point(
2152+
sys, 0, 0, root_method='lm', root_kwargs={'tol': 1e-6})
2153+
2154+
# Make sure we can get back the right arguments in a tuple
2155+
op_point = ct.find_operating_point(sys, 0, 0, return_outputs=True)
2156+
assert len(op_point) == 3
2157+
assert isinstance(op_point[0], np.ndarray)
2158+
assert isinstance(op_point[1], np.ndarray)
2159+
assert isinstance(op_point[2], np.ndarray)
2160+
2161+
with pytest.warns(FutureWarning, match="return_outputs"):
2162+
op_point = ct.find_operating_point(sys, 0, 0, return_y=True)
2163+
assert len(op_point) == 3
2164+
assert isinstance(op_point[0], np.ndarray)
2165+
assert isinstance(op_point[1], np.ndarray)
2166+
assert isinstance(op_point[2], np.ndarray)
2167+
2168+
# Make sure we can get back the right arguments in a tuple
2169+
op_point = ct.find_operating_point(sys, 0, 0, return_result=True)
2170+
assert len(op_point) == 3
2171+
assert isinstance(op_point[0], np.ndarray)
2172+
assert isinstance(op_point[1], np.ndarray)
2173+
assert isinstance(op_point[2], scipy.optimize.OptimizeResult)
2174+
2175+
# Make sure we can get back the right arguments in a tuple
2176+
op_point = ct.find_operating_point(
2177+
sys, 0, 0, return_result=True, return_outputs=True)
2178+
assert len(op_point) == 4
2179+
assert isinstance(op_point[0], np.ndarray)
2180+
assert isinstance(op_point[1], np.ndarray)
2181+
assert isinstance(op_point[2], np.ndarray)
2182+
assert isinstance(op_point[3], scipy.optimize.OptimizeResult)
2183+
2184+
20902185
def test_iosys_sample():
20912186
csys = ct.rss(2, 1, 1)
20922187
dsys = csys.sample(0.1)

‎control/tests/kwargs_test.py

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@
2424
import control.tests.frd_test as frd_test
2525
import control.tests.freqplot_test as freqplot_test
2626
import control.tests.interconnect_test as interconnect_test
27+
import control.tests.iosys_test as iosys_test
2728
import control.tests.optimal_test as optimal_test
2829
import control.tests.statefbk_test as statefbk_test
2930
import control.tests.stochsys_test as stochsys_test
@@ -252,6 +253,8 @@ def test_response_plot_kwargs(data_fcn, plot_fcn, mimo):
252253
'dlqr': test_unrecognized_kwargs,
253254
'drss': test_unrecognized_kwargs,
254255
'feedback': test_unrecognized_kwargs,
256+
'find_eqpt': iosys_test.test_find_operating_point,
257+
'find_operating_point': iosys_test.test_find_operating_point,
255258
'flatsys.flatsys': test_unrecognized_kwargs,
256259
'frd': frd_test.TestFRD.test_unrecognized_keyword,
257260
'gangof4': test_matplotlib_kwargs,

‎doc/control.rst

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -144,7 +144,7 @@ Nonlinear system support
144144
:toctree: generated/
145145

146146
describing_function
147-
find_eqpt
147+
find_operating_point
148148
linearize
149149
input_output_response
150150
summing_junction

‎doc/conventions.rst

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -289,10 +289,10 @@ element of the `control.config.defaults` dictionary::
289289

290290
ct.config.defaults['module.parameter'] = value
291291

292-
The `~control.config.set_defaults` function can also be used to set multiple
293-
configuration parameters at the same time::
292+
The :func:`~control.set_defaults` function can also be used to
293+
set multiple configuration parameters at the same time::
294294

295-
ct.config.set_defaults('module', param1=val1, param2=val2, ...]
295+
ct.set_defaults('module', param1=val1, param2=val2, ...]
296296

297297
Finally, there are also functions available set collections of variables based
298298
on standard configurations.

‎doc/intro.rst

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -56,7 +56,7 @@ they are not already present.
5656
.. note::
5757
Mixing packages from conda-forge and the default conda channel
5858
can sometimes cause problems with dependencies, so it is usually best to
59-
instally NumPy, SciPy, and Matplotlib from conda-forge as well.
59+
install NumPy, SciPy, and Matplotlib from conda-forge as well.
6060

6161
To install using pip::
6262

‎doc/iosys.rst

Lines changed: 10 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -16,12 +16,13 @@ function::
1616
resp = ct.input_output_response(io_sys, T, U, X0, params)
1717
t, y, x = resp.time, resp.outputs, resp.states
1818

19-
An input/output system can be linearized around an equilibrium point to obtain
20-
a :class:`~control.StateSpace` linear system. Use the
21-
:func:`~control.find_eqpt` function to obtain an equilibrium point and the
22-
:func:`~control.linearize` function to linearize about that equilibrium point::
19+
An input/output system can be linearized around an equilibrium point
20+
to obtain a :class:`~control.StateSpace` linear system. Use the
21+
:func:`~control.find_operating_point` function to obtain an
22+
equilibrium point and the :func:`~control.linearize` function to
23+
linearize about that equilibrium point::
2324

24-
xeq, ueq = ct.find_eqpt(io_sys, X0, U0)
25+
xeq, ueq = ct.find_operating_point(io_sys, X0, U0)
2526
ss_sys = ct.linearize(io_sys, xeq, ueq)
2627

2728
Input/output systems are automatically created for state space LTI systems
@@ -123,9 +124,8 @@ system and computing the linearization about that point.
123124

124125
.. code-block:: python
125126
126-
eqpt = ct.find_eqpt(io_predprey, X0, 0)
127-
xeq = eqpt[0] # choose the nonzero equilibrium point
128-
lin_predprey = ct.linearize(io_predprey, xeq, 0)
127+
eqpt = ct.find_operating_point(io_predprey, X0, 0)
128+
lin_predprey = ct.linearize(io_predprey, eqpt)
129129
130130
We next compute a controller that stabilizes the equilibrium point using
131131
eigenvalue placement and computing the feedforward gain using the number of
@@ -544,11 +544,12 @@ Module classes and functions
544544
~control.InterconnectedSystem
545545
~control.LinearICSystem
546546
~control.NonlinearIOSystem
547+
~control.OperatingPoint
547548

548549
.. autosummary::
549550
:toctree: generated/
550551

551-
~control.find_eqpt
552+
~control.find_operating_point
552553
~control.interconnect
553554
~control.input_output_response
554555
~control.linearize

‎examples/cruise-control.py

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -165,7 +165,7 @@ def motor_torque(omega, params={}):
165165

166166
for m in (1200, 1600, 2000):
167167
# Compute the equilibrium state for the system
168-
X0, U0 = ct.find_eqpt(
168+
X0, U0 = ct.find_operating_point(
169169
cruise_tf, [0, vref[0]], [vref[0], gear[0], theta0[0]],
170170
iu=[1, 2], y0=[vref[0], 0], iy=[0], params={'m': m})
171171

@@ -347,9 +347,9 @@ def cruise_plot(sys, t, y, label=None, t_hill=None, vref=20, antiwindup=False,
347347

348348
# Compute the equilibrium throttle setting for the desired speed (solve for x
349349
# and u given the gear, slope, and desired output velocity)
350-
X0, U0, Y0 = ct.find_eqpt(
350+
X0, U0, Y0 = ct.find_operating_point(
351351
cruise_pi, [vref[0], 0], [vref[0], gear[0], theta0[0]],
352-
y0=[0, vref[0]], iu=[1, 2], iy=[1], return_y=True)
352+
y0=[0, vref[0]], iu=[1, 2], iy=[1], return_outputs=True)
353353

354354
# Now simulate the effect of a hill at t = 5 seconds
355355
plt.figure()

‎examples/cruise.ipynb

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -804,7 +804,7 @@
804804
"# Compute the equilibrium throttle setting for the desired speed\n",
805805
"X0, U0, Y0 = ct.find_eqpt(\n",
806806
" cruise_pi, [vref[0], 0], [vref[0], gear[0], theta0[0]],\n",
807-
" y0=[0, vref[0]], iu=[1, 2], iy=[1], return_y=True)\n",
807+
" y0=[0, vref[0]], iu=[1, 2], iy=[1], return_outputs=True)\n",
808808
"\n",
809809
"# Now simulate the effect of a hill at t = 5 seconds\n",
810810
"plt.figure()\n",

0 commit comments

Comments
 (0)
0