Skip to content

Commit 71bd731

Browse files
murrayrmmurrayrm
authored
Merge pull request #1081 from sdahdah/main
Support for binary operations between MIMO and SISO LTI systems
2 parents bddc796 + ccf9ce1 commit 71bd731

File tree

10 files changed

+1138
-86
lines changed

10 files changed

+1138
-86
lines changed

‎control/bdalg.py‎

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -356,14 +356,14 @@ def feedback(sys1, sys2=1, sign=-1, **kwargs):
356356
def append(*sys, **kwargs):
357357
"""append(sys1, sys2, [..., sysn])
358358
359-
Group LTI state space models by appending their inputs and outputs.
359+
Group LTI models by appending their inputs and outputs.
360360
361361
Forms an augmented system model, and appends the inputs and
362362
outputs together.
363363
364364
Parameters
365365
----------
366-
sys1, sys2, ..., sysn: scalar, array, or :class:`StateSpace`
366+
sys1, sys2, ..., sysn : scalar, array, or :class:`LTI`
367367
I/O systems to combine.
368368
369369
Other Parameters
@@ -382,9 +382,10 @@ def append(*sys, **kwargs):
382382
383383
Returns
384384
-------
385-
out: :class:`StateSpace`
385+
out : :class:`LTI`
386386
Combined system, with input/output vectors consisting of all
387-
input/output vectors appended.
387+
input/output vectors appended. Specific type returned is the type of
388+
the first argument.
388389
389390
See Also
390391
--------
@@ -405,7 +406,7 @@ def append(*sys, **kwargs):
405406
(3, 8, 7)
406407
407408
"""
408-
s1 = ss._convert_to_statespace(sys[0])
409+
s1 = sys[0]
409410
for s in sys[1:]:
410411
s1 = s1.append(s)
411412
s1.update_names(**kwargs)

‎control/frdata.py‎

Lines changed: 45 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@
2020
from scipy.interpolate import splev, splprep
2121

2222
from . import config
23+
from . import bdalg
2324
from .exception import pandas_check
2425
from .iosys import InputOutputSystem, NamedSignal, _extended_system_name, \
2526
_process_iosys_keywords, _process_subsys_index, common_timebase
@@ -455,6 +456,12 @@ def __add__(self, other):
455456
else:
456457
other = _convert_to_frd(other, omega=self.omega)
457458

459+
# Promote SISO object to compatible dimension
460+
if self.issiso() and not other.issiso():
461+
self = np.ones((other.noutputs, other.ninputs)) * self
462+
elif not self.issiso() and other.issiso():
463+
other = np.ones((self.noutputs, self.ninputs)) * other
464+
458465
# Check that the input-output sizes are consistent.
459466
if self.ninputs != other.ninputs:
460467
raise ValueError(
@@ -492,6 +499,12 @@ def __mul__(self, other):
492499
else:
493500
other = _convert_to_frd(other, omega=self.omega)
494501

502+
# Promote SISO object to compatible dimension
503+
if self.issiso() and not other.issiso():
504+
self = bdalg.append(*([self] * other.noutputs))
505+
elif not self.issiso() and other.issiso():
506+
other = bdalg.append(*([other] * self.ninputs))
507+
495508
# Check that the input-output sizes are consistent.
496509
if self.ninputs != other.noutputs:
497510
raise ValueError(
@@ -519,6 +532,12 @@ def __rmul__(self, other):
519532
else:
520533
other = _convert_to_frd(other, omega=self.omega)
521534

535+
# Promote SISO object to compatible dimension
536+
if self.issiso() and not other.issiso():
537+
self = bdalg.append(*([self] * other.ninputs))
538+
elif not self.issiso() and other.issiso():
539+
other = bdalg.append(*([other] * self.noutputs))
540+
522541
# Check that the input-output sizes are consistent.
523542
if self.noutputs != other.ninputs:
524543
raise ValueError(
@@ -547,11 +566,9 @@ def __truediv__(self, other):
547566
else:
548567
other = _convert_to_frd(other, omega=self.omega)
549568

550-
if (self.ninputs > 1 or self.noutputs > 1 or
551-
other.ninputs > 1 or other.noutputs > 1):
552-
raise NotImplementedError(
553-
"FRD.__truediv__ is currently only implemented for SISO "
554-
"systems.")
569+
if (other.ninputs > 1 or other.noutputs > 1):
570+
# FRD.__truediv__ is currently only implemented for SISO systems
571+
return NotImplemented
555572

556573
return FRD(self.fresp/other.fresp, self.omega,
557574
smooth=(self._ifunc is not None) and
@@ -566,11 +583,9 @@ def __rtruediv__(self, other):
566583
else:
567584
other = _convert_to_frd(other, omega=self.omega)
568585

569-
if (self.ninputs > 1 or self.noutputs > 1 or
570-
other.ninputs > 1 or other.noutputs > 1):
571-
raise NotImplementedError(
572-
"FRD.__rtruediv__ is currently only implemented for "
573-
"SISO systems.")
586+
if (self.ninputs > 1 or self.noutputs > 1):
587+
# FRD.__rtruediv__ is currently only implemented for SISO systems
588+
return NotImplemented
574589

575590
return other / self
576591

@@ -803,6 +818,26 @@ def feedback(self, other=1, sign=-1):
803818

804819
return FRD(fresp, other.omega, smooth=(self._ifunc is not None))
805820

821+
def append(self, other):
822+
"""Append a second model to the present model.
823+
824+
The second model is converted to FRD if necessary, inputs and
825+
outputs are appended and their order is preserved"""
826+
other = _convert_to_frd(other, omega=self.omega, inputs=other.ninputs,
827+
outputs=other.noutputs)
828+
829+
# TODO: handle omega re-mapping
830+
831+
new_fresp = np.zeros(
832+
(self.noutputs + other.noutputs, self.ninputs + other.ninputs,
833+
self.omega.shape[-1]), dtype=complex)
834+
new_fresp[:self.noutputs, :self.ninputs, :] = np.reshape(
835+
self.fresp, (self.noutputs, self.ninputs, -1))
836+
new_fresp[self.noutputs:, self.ninputs:, :] = np.reshape(
837+
other.fresp, (other.noutputs, other.ninputs, -1))
838+
839+
return FRD(new_fresp, self.omega, smooth=(self._ifunc is not None))
840+
806841
# Plotting interface
807842
def plot(self, plot_type=None, *args, **kwargs):
808843
"""Plot the frequency response using a Bode plot.

‎control/statesp.py‎

Lines changed: 70 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@
3030
from scipy.signal import cont2discrete
3131

3232
from . import config
33+
from . import bdalg
3334
from .exception import ControlMIMONotImplemented, ControlSlycot, slycot_check
3435
from .frdata import FrequencyResponseData
3536
from .iosys import InputOutputSystem, NamedSignal, _process_dt_keyword, \
@@ -572,6 +573,9 @@ def __add__(self, other):
572573

573574
elif isinstance(other, np.ndarray):
574575
other = np.atleast_2d(other)
576+
# Special case for SISO
577+
if self.issiso():
578+
self = np.ones_like(other) * self
575579
if self.ninputs != other.shape[0]:
576580
raise ValueError("array has incompatible shape")
577581
A, B, C = self.A, self.B, self.C
@@ -582,6 +586,12 @@ def __add__(self, other):
582586
return NotImplemented # let other.__rmul__ handle it
583587

584588
else:
589+
# Promote SISO object to compatible dimension
590+
if self.issiso() and not other.issiso():
591+
self = np.ones((other.noutputs, other.ninputs)) * self
592+
elif not self.issiso() and other.issiso():
593+
other = np.ones((self.noutputs, self.ninputs)) * other
594+
585595
# Check to make sure the dimensions are OK
586596
if ((self.ninputs != other.ninputs) or
587597
(self.noutputs != other.noutputs)):
@@ -636,6 +646,10 @@ def __mul__(self, other):
636646

637647
elif isinstance(other, np.ndarray):
638648
other = np.atleast_2d(other)
649+
# Special case for SISO
650+
if self.issiso():
651+
self = bdalg.append(*([self] * other.shape[0]))
652+
# Dimension check after broadcasting
639653
if self.ninputs != other.shape[0]:
640654
raise ValueError("array has incompatible shape")
641655
A, C = self.A, self.C
@@ -647,6 +661,12 @@ def __mul__(self, other):
647661
return NotImplemented # let other.__rmul__ handle it
648662

649663
else:
664+
# Promote SISO object to compatible dimension
665+
if self.issiso() and not other.issiso():
666+
self = bdalg.append(*([self] * other.noutputs))
667+
elif not self.issiso() and other.issiso():
668+
other = bdalg.append(*([other] * self.ninputs))
669+
650670
# Check to make sure the dimensions are OK
651671
if self.ninputs != other.noutputs:
652672
raise ValueError(
@@ -686,23 +706,67 @@ def __rmul__(self, other):
686706
return StateSpace(self.A, B, self.C, D, self.dt)
687707

688708
elif isinstance(other, np.ndarray):
689-
C = np.atleast_2d(other) @ self.C
690-
D = np.atleast_2d(other) @ self.D
709+
other = np.atleast_2d(other)
710+
# Special case for SISO transfer function
711+
if self.issiso():
712+
self = bdalg.append(*([self] * other.shape[1]))
713+
# Dimension check after broadcasting
714+
if self.noutputs != other.shape[1]:
715+
raise ValueError("array has incompatible shape")
716+
C = other @ self.C
717+
D = other @ self.D
691718
return StateSpace(self.A, self.B, C, D, self.dt)
692719

693720
if not isinstance(other, StateSpace):
694721
return NotImplemented
695722

723+
# Promote SISO object to compatible dimension
724+
if self.issiso() and not other.issiso():
725+
self = bdalg.append(*([self] * other.ninputs))
726+
elif not self.issiso() and other.issiso():
727+
other = bdalg.append(*([other] * self.noutputs))
728+
696729
return other * self
697730

698731
# TODO: general __truediv__ requires descriptor system support
699732
def __truediv__(self, other):
700733
"""Division of state space systems by TFs, FRDs, scalars, and arrays"""
701-
if not isinstance(other, (LTI, InputOutputSystem)):
702-
return self * (1/other)
703-
else:
734+
# Let ``other.__rtruediv__`` handle it
735+
try:
736+
return self * (1 / other)
737+
except ValueError:
704738
return NotImplemented
705739

740+
def __rtruediv__(self, other):
741+
"""Division by state space system"""
742+
return other * self**-1
743+
744+
def __pow__(self, other):
745+
"""Power of a state space system"""
746+
if not type(other) == int:
747+
raise ValueError("Exponent must be an integer")
748+
if self.ninputs != self.noutputs:
749+
# System must have same number of inputs and outputs
750+
return NotImplemented
751+
if other < -1:
752+
return (self**-1)**(-other)
753+
elif other == -1:
754+
try:
755+
Di = scipy.linalg.inv(self.D)
756+
except scipy.linalg.LinAlgError:
757+
# D matrix must be nonsingular
758+
return NotImplemented
759+
Ai = self.A - self.B @ Di @ self.C
760+
Bi = self.B @ Di
761+
Ci = -Di @ self.C
762+
return StateSpace(Ai, Bi, Ci, Di, self.dt)
763+
elif other == 0:
764+
return StateSpace([], [], [], np.eye(self.ninputs), self.dt)
765+
elif other == 1:
766+
return self
767+
elif other > 1:
768+
return self * (self**(other - 1))
769+
706770
def __call__(self, x, squeeze=None, warn_infinite=True):
707771
"""Evaluate system's frequency response at complex frequencies.
708772
@@ -1107,7 +1171,7 @@ def minreal(self, tol=0.0):
11071171
A, B, C, nr = tb01pd(self.nstates, self.ninputs, self.noutputs,
11081172
self.A, B, C, tol=tol)
11091173
return StateSpace(A[:nr, :nr], B[:nr, :self.ninputs],
1110-
C[:self.noutputs, :nr], self.D)
1174+
C[:self.noutputs, :nr], self.D, self.dt)
11111175
except ImportError:
11121176
raise TypeError("minreal requires slycot tb01pd")
11131177
else:

‎control/tests/bdalg_test.py‎

Lines changed: 1 addition & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@
88
import pytest
99

1010
import control as ctrl
11-
from control.xferfcn import TransferFunction
11+
from control.xferfcn import TransferFunction, _tf_close_coeff
1212
from control.statesp import StateSpace
1313
from control.bdalg import feedback, append, connect
1414
from control.lti import zeros, poles
@@ -870,50 +870,3 @@ def test_error_combine_tf(self, tf_array, exception):
870870
"""Test error cases."""
871871
with pytest.raises(exception):
872872
ctrl.combine_tf(tf_array)
873-
874-
875-
def _tf_close_coeff(tf_a, tf_b, rtol=1e-5, atol=1e-8):
876-
"""Check if two transfer functions have close coefficients.
877-
878-
Parameters
879-
----------
880-
tf_a : TransferFunction
881-
First transfer function.
882-
tf_b : TransferFunction
883-
Second transfer function.
884-
rtol : float
885-
Relative tolerance for ``np.allclose``.
886-
atol : float
887-
Absolute tolerance for ``np.allclose``.
888-
889-
Returns
890-
-------
891-
bool
892-
True if transfer function cofficients are all close.
893-
"""
894-
# Check number of outputs and inputs
895-
if tf_a.noutputs != tf_b.noutputs:
896-
return False
897-
if tf_a.ninputs != tf_b.ninputs:
898-
return False
899-
# Check timestep
900-
if tf_a.dt != tf_b.dt:
901-
return False
902-
# Check coefficient arrays
903-
for i in range(tf_a.noutputs):
904-
for j in range(tf_a.ninputs):
905-
if not np.allclose(
906-
tf_a.num_array[i, j],
907-
tf_b.num_array[i, j],
908-
rtol=rtol,
909-
atol=atol,
910-
):
911-
return False
912-
if not np.allclose(
913-
tf_a.den_array[i, j],
914-
tf_b.den_array[i, j],
915-
rtol=rtol,
916-
atol=atol,
917-
):
918-
return False
919-
return True

‎control/tests/docstrings_test.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,7 @@
2828

2929
# Checksums to use for checking whether a docstring has changed
3030
function_docstring_hash = {
31-
control.append: '48548c4c4e0083312b3ea9e56174b0b5',
31+
control.append: '1bddbac0fe932755c85e9fb0bfb97d88',
3232
control.describing_function_plot: '95f894706b1d3eeb3b854934596af09f',
3333
control.dlqe: '9db995ed95c2214ce97074b0616a3191',
3434
control.dlqr: '896cfa651dbbd80e417635904d13c9d6',
@@ -37,7 +37,7 @@
3737
control.margin: 'f02b3034f5f1d44ce26f916cc3e51600',
3838
control.parallel: 'bfc470aef75dbb923f9c6fb8bf3c9b43',
3939
control.series: '9aede1459667738f05cf4fc46603a4f6',
40-
control.ss2tf: '48ff25d22d28e7b396e686dd5eb58831',
40+
control.ss2tf: 'e779b8d70205bc1218cc2a4556a66e4b',
4141
control.tf2ss: '086a3692659b7321c2af126f79f4bc11',
4242
control.markov: 'a4199c54cb50f07c0163d3790739eafe',
4343
control.gangof4: '0e52eb6cf7ce024f9a41f3ae3ebf04f7',

0 commit comments

Comments
 (0)