Skip to content

fixed prewarp not working in c2d and sample_system, margin docstring improvements #669

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 4 commits into from
Nov 8, 2021
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 18 additions & 13 deletions control/dtime.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
Routines in this module:

sample_system()
c2d()
"""

"""Copyright (c) 2012 by California Institute of Technology
Expand Down Expand Up @@ -58,16 +59,19 @@ def sample_system(sysc, Ts, method='zoh', alpha=None, prewarp_frequency=None):

Parameters
----------
sysc : LTI (StateSpace or TransferFunction)
sysc : LTI (:class:`StateSpace` or :class:`TransferFunction`)
Continuous time system to be converted
Ts : real > 0
Ts : float > 0
Sampling period
method : string
Method to use for conversion, e.g. 'bilinear', 'zoh' (default)

prewarp_frequency : real within [0, infinity)
alpha : float within [0, 1]
The generalized bilinear transformation weighting parameter, which
should only be specified with method="gbt", and is ignored
otherwise. See :func:`scipy.signal.cont2discrete`.
prewarp_frequency : float within [0, infinity)
The frequency [rad/s] at which to match with the input continuous-
time system's magnitude and phase
time system's magnitude and phase (only valid for method='bilinear')

Returns
-------
Expand All @@ -76,7 +80,7 @@ def sample_system(sysc, Ts, method='zoh', alpha=None, prewarp_frequency=None):

Notes
-----
See :meth:`StateSpace.sample` or :meth:`TransferFunction.sample`` for
See :meth:`StateSpace.sample` or :meth:`TransferFunction.sample` for
further details.

Examples
Expand All @@ -89,7 +93,8 @@ def sample_system(sysc, Ts, method='zoh', alpha=None, prewarp_frequency=None):
if not isctime(sysc):
raise ValueError("First argument must be continuous time system")

return sysc.sample(Ts, method, alpha, prewarp_frequency)
return sysc.sample(Ts,
method=method, alpha=alpha, prewarp_frequency=prewarp_frequency)


def c2d(sysc, Ts, method='zoh', prewarp_frequency=None):
Expand All @@ -98,20 +103,19 @@ def c2d(sysc, Ts, method='zoh', prewarp_frequency=None):

Parameters
----------
sysc : LTI (StateSpace or TransferFunction)
sysc : LTI (:class:`StateSpace` or :class:`TransferFunction`)
Continuous time system to be converted
Ts : real > 0
Ts : float > 0
Sampling period
method : string
Method to use for conversion, e.g. 'bilinear', 'zoh' (default)

prewarp_frequency : real within [0, infinity)
The frequency [rad/s] at which to match with the input continuous-
time system's magnitude and phase
time system's magnitude and phase (only valid for method='bilinear')

Returns
-------
sysd : linsys
sysd : LTI of the same class
Discrete time system, with sampling rate Ts

Notes
Expand All @@ -126,6 +130,7 @@ def c2d(sysc, Ts, method='zoh', prewarp_frequency=None):
"""

# Call the sample_system() function to do the work
sysd = sample_system(sysc, Ts, method, prewarp_frequency)
sysd = sample_system(sysc, Ts,
method=method, prewarp_frequency=prewarp_frequency)

return sysd
20 changes: 12 additions & 8 deletions control/margins.py
Original file line number Diff line number Diff line change
Expand Up @@ -283,14 +283,16 @@ def stability_margins(sysdata, returnall=False, epsw=0.0, method='best'):
-------
gm : float or array_like
Gain margin
pm : float or array_loke
pm : float or array_like
Phase margin
sm : float or array_like
Stability margin, the minimum distance from the Nyquist plot to -1
wpc : float or array_like
Phase crossover frequency (where phase crosses -180 degrees)
Phase crossover frequency (where phase crosses -180 degrees), which is
associated with the gain margin.
wgc : float or array_like
Gain crossover frequency (where gain crosses 1)
Gain crossover frequency (where gain crosses 1), which is associated
with the phase margin.
wms : float or array_like
Stability margin frequency (where Nyquist plot is closest to -1)

Expand Down Expand Up @@ -522,10 +524,12 @@ def margin(*args):
Gain margin
pm : float
Phase margin (in degrees)
wpc : float or array_like
Phase crossover frequency (where phase crosses -180 degrees)
wgc : float or array_like
Gain crossover frequency (where gain crosses 1)
wcg : float or array_like
Crossover frequency associated with gain margin (phase crossover
frequency), where phase crosses below -180 degrees.
wcp : float or array_like
Crossover frequency associated with phase margin (gain crossover
frequency), where gain crosses below 1.

Margins are calculated for a SISO open-loop system.

Expand All @@ -536,7 +540,7 @@ def margin(*args):
Examples
--------
>>> sys = tf(1, [1, 2, 1, 0])
>>> gm, pm, wg, wp = margin(sys)
>>> gm, pm, wcg, wcp = margin(sys)

"""
if len(args) == 1:
Expand Down
18 changes: 14 additions & 4 deletions control/tests/discrete_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,8 @@
import pytest

from control import (StateSpace, TransferFunction, bode, common_timebase,
evalfr, feedback, forced_response, impulse_response,
isctime, isdtime, rss, sample_system, step_response,
feedback, forced_response, impulse_response,
isctime, isdtime, rss, c2d, sample_system, step_response,
timebase)


Expand Down Expand Up @@ -382,10 +382,20 @@ def test_sample_system_prewarp(self, tsys, plantname):
Ts = 0.025
# test state space version
plant = getattr(tsys, plantname)
plant_fr = plant(wwarp * 1j)

plant_d_warped = plant.sample(Ts, 'bilinear', prewarp_frequency=wwarp)
plant_fr = evalfr(plant, wwarp * 1j)
dt = plant_d_warped.dt
plant_d_fr = evalfr(plant_d_warped, np.exp(wwarp * 1.j * dt))
plant_d_fr = plant_d_warped(np.exp(wwarp * 1.j * dt))
np.testing.assert_array_almost_equal(plant_fr, plant_d_fr)

plant_d_warped = sample_system(plant, Ts, 'bilinear',
prewarp_frequency=wwarp)
plant_d_fr = plant_d_warped(np.exp(wwarp * 1.j * dt))
np.testing.assert_array_almost_equal(plant_fr, plant_d_fr)

plant_d_warped = c2d(plant, Ts, 'bilinear', prewarp_frequency=wwarp)
plant_d_fr = plant_d_warped(np.exp(wwarp * 1.j * dt))
np.testing.assert_array_almost_equal(plant_fr, plant_d_fr)

def test_sample_system_errors(self, tsys):
Expand Down
20 changes: 10 additions & 10 deletions control/tests/matlab_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -355,13 +355,13 @@ def testLsim_mimo(self, mimo):
def testMargin(self, siso):
"""Test margin()"""
#! TODO: check results to make sure they are OK
gm, pm, wg, wp = margin(siso.tf1)
gm, pm, wg, wp = margin(siso.tf2)
gm, pm, wg, wp = margin(siso.ss1)
gm, pm, wg, wp = margin(siso.ss2)
gm, pm, wg, wp = margin(siso.ss2 * siso.ss2 * 2)
gm, pm, wcg, wcp = margin(siso.tf1)
gm, pm, wcg, wcp = margin(siso.tf2)
gm, pm, wcg, wcp = margin(siso.ss1)
gm, pm, wcg, wcp = margin(siso.ss2)
gm, pm, wcg, wcp = margin(siso.ss2 * siso.ss2 * 2)
np.testing.assert_array_almost_equal(
[gm, pm, wg, wp], [1.5451, 75.9933, 1.2720, 0.6559], decimal=3)
[gm, pm, wcg, wcp], [1.5451, 75.9933, 1.2720, 0.6559], decimal=3)

def testDcgain(self, siso):
"""Test dcgain() for SISO system"""
Expand Down Expand Up @@ -781,12 +781,12 @@ def testCombi01(self):
# total open loop
Hol = Hc*Hno*Hp

gm, pm, wg, wp = margin(Hol)
# print("%f %f %f %f" % (gm, pm, wg, wp))
gm, pm, wcg, wcp = margin(Hol)
# print("%f %f %f %f" % (gm, pm, wcg, wcp))
np.testing.assert_allclose(gm, 3.32065569155)
np.testing.assert_allclose(pm, 46.9740430224)
np.testing.assert_allclose(wg, 0.176469728448)
np.testing.assert_allclose(wp, 0.0616288455466)
np.testing.assert_allclose(wcg, 0.176469728448)
np.testing.assert_allclose(wcp, 0.0616288455466)

def test_tf_string_args(self):
"""Make sure s and z are defined properly"""
Expand Down
18 changes: 9 additions & 9 deletions control/xferfcn.py
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,7 @@
from itertools import chain
from re import sub
from .lti import LTI, common_timebase, isdtime, _process_frequency_response
from .exception import ControlMIMONotImplemented
from . import config

__all__ = ['TransferFunction', 'tf', 'ss2tf', 'tfdata']
Expand Down Expand Up @@ -793,9 +794,9 @@ def feedback(self, other=1, sign=-1):
if (self.ninputs > 1 or self.noutputs > 1 or
other.ninputs > 1 or other.noutputs > 1):
# TODO: MIMO feedback
raise NotImplementedError(
"TransferFunction.feedback is currently only implemented "
"for SISO functions.")
raise ControlMIMONotImplemented(
"TransferFunction.feedback is currently not implemented for "
"MIMO systems.")
dt = common_timebase(self.dt, other.dt)

num1 = self.num[0][0]
Expand Down Expand Up @@ -1085,12 +1086,10 @@ def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None):
* euler: Euler (or forward difference) method ("gbt" with alpha=0)
* backward_diff: Backwards difference ("gbt" with alpha=1.0)
* zoh: zero-order hold (default)

alpha : float within [0, 1]
The generalized bilinear transformation weighting parameter, which
should only be specified with method="gbt", and is ignored
otherwise.

otherwise. See :func:`scipy.signal.cont2discrete`.
prewarp_frequency : float within [0, infinity)
The frequency [rad/s] at which to match with the input continuous-
time system's magnitude and phase (the gain=1 crossover frequency,
Expand All @@ -1100,7 +1099,7 @@ def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None):
Returns
-------
sysd : TransferFunction system
Discrete time system, with sampling rate Ts
Discrete time system, with sample period Ts

Notes
-----
Expand All @@ -1117,7 +1116,7 @@ def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None):
if not self.isctime():
raise ValueError("System must be continuous time system")
if not self.issiso():
raise NotImplementedError("MIMO implementation not available")
raise ControlMIMONotImplemented("Not implemented for MIMO systems")
if method == "matched":
return _c2d_matched(self, Ts)
sys = (self.num[0][0], self.den[0][0])
Expand Down Expand Up @@ -1373,7 +1372,8 @@ def _convert_to_transfer_function(sys, **kw):
except ImportError:
# If slycot is not available, use signal.lti (SISO only)
if sys.ninputs != 1 or sys.noutputs != 1:
raise TypeError("No support for MIMO without slycot.")
raise ControlMIMONotImplemented("Not implemented for " +
"MIMO systems without slycot.")

# Do the conversion using sp.signal.ss2tf
# Note that this returns a 2D array for the numerator
Expand Down