Repository navigation
Expand file tree
/
Copy pathnlsys.py
More file actions
2808 lines (2333 loc) · 115 KB
/
Copy pathnlsys.py
File metadata and controls
2808 lines (2333 loc) · 115 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# nlsys.py - input/output system module
# RMM, 28 April 2019
#
# Additional features to add
# * Allow constant inputs for MIMO input_output_response (w/out ones)
# * Add support for constants/matrices as part of operators (1 + P)
# * Add unit tests (and example?) for time-varying systems
# * Allow time vector for discrete time simulations to be multiples of dt
# * Check the way initial outputs for discrete time systems are handled
#
"""The :mod:`~control.nlsys` module contains the
:class:`~control.NonlinearIOSystem` class that represents (possibly nonlinear)
input/output systems. The :class:`~control.NonlinearIOSystem` class is a
general class that defines any continuous or discrete time dynamical system.
Input/output systems can be simulated and also used to compute operating
points and linearizations.
"""
import copy
from warnings import warn
import numpy as np
import scipy as sp
from . import config
from .iosys import InputOutputSystem, _parse_spec, _process_iosys_keywords, \
_process_signal_list, common_timebase, isctime, isdtime
from .timeresp import _check_convert_array, _process_time_response, \
TimeResponseData, TimeResponseList
__all__ = ['NonlinearIOSystem', 'InterconnectedSystem', 'nlsys',
'input_output_response', 'find_eqpt', 'linearize',
'interconnect', 'connection_table', 'OperatingPoint',
'find_operating_point']
class NonlinearIOSystem(InputOutputSystem):
"""Nonlinear I/O system.
Creates an :class:`~control.InputOutputSystem` for a nonlinear system by
specifying a state update function and an output function. The new system
can be a continuous or discrete time system (Note: discrete-time systems
are not yet supported by most functions.)
Parameters
----------
updfcn : callable
Function returning the state update function
`updfcn(t, x, u, params) -> array`
where `x` is a 1-D array with shape (nstates,), `u` is a 1-D array
with shape (ninputs,), `t` is a float representing the currrent
time, and `params` is a dict containing the values of parameters
used by the function.
outfcn : callable
Function returning the output at the given state
`outfcn(t, x, u, params) -> array`
where the arguments are the same as for `upfcn`.
inputs : int, list of str or None, optional
Description of the system inputs. This can be given as an integer
count or as a list of strings that name the individual signals.
If an integer count is specified, the names of the signal will be
of the form 's[i]' (where 's' is one of 'u', 'y', or 'x'). If
this parameter is not given or given as `None`, the relevant
quantity will be determined when possible based on other
information provided to functions using the system.
outputs : int, list of str or None, optional
Description of the system outputs. Same format as `inputs`.
states : int, list of str, or None, optional
Description of the system states. Same format as `inputs`.
dt : timebase, optional
The timebase for the system, used to specify whether the system is
operating in continuous or discrete time. It can have the
following values:
* dt = 0: continuous time system (default)
* dt > 0: discrete time system with sampling period 'dt'
* dt = True: discrete time with unspecified sampling period
* dt = None: no timebase specified
name : string, optional
System name (used for specifying signals). If unspecified, a
generic name <sys[id]> is generated with a unique integer id.
params : dict, optional
Parameter values for the system. Passed to the evaluation functions
for the system as default values, overriding internal defaults.
See Also
--------
InputOutputSystem : Input/output system class.
Notes
-----
The :class:`~control.InputOuputSystem` class (and its subclasses) makes
use of two special methods for implementing much of the work of the class:
* _rhs(t, x, u): compute the right hand side of the differential or
difference equation for the system. If not specified, the system
has no state.
* _out(t, x, u): compute the output for the current state of the system.
The default is to return the entire system state.
"""
def __init__(self, updfcn, outfcn=None, params=None, **kwargs):
"""Create a nonlinear I/O system given update and output functions."""
# Process keyword arguments
name, inputs, outputs, states, dt = _process_iosys_keywords(kwargs)
# Initialize the rest of the structure
super().__init__(
inputs=inputs, outputs=outputs, states=states, dt=dt, name=name,
**kwargs
)
self.params = {} if params is None else params.copy()
# Store the update and output functions
self.updfcn = updfcn
self.outfcn = outfcn
# Check to make sure arguments are consistent
if updfcn is None:
if self.nstates is None:
self.nstates = 0
self.updfcn = lambda t, x, u, params: np.zeros(0)
else:
raise ValueError(
"states specified but no update function given.")
if outfcn is None:
if self.noutputs == 0:
self.outfcn = lambda t, x, u, params: np.zeros(0)
elif self.noutputs is None and self.nstates is not None:
self.noutputs = self.nstates
if len(self.output_index) == 0:
# Use state names for outputs
self.output_index = self.state_index
elif self.noutputs is not None and self.noutputs == self.nstates:
# Number of outputs = number of states => all is OK
pass
elif self.noutputs is not None and self.noutputs != 0:
raise ValueError("outputs specified but no output function "
"(and nstates not known).")
# Initialize current parameters to default parameters
self._current_params = {} if params is None else params.copy()
def __str__(self):
return f"{InputOutputSystem.__str__(self)}\n\n" + \
f"Update: {self.updfcn}\n" + \
f"Output: {self.outfcn}"
# Return the value of a static nonlinear system
def __call__(sys, u, params=None, squeeze=None):
"""Evaluate a (static) nonlinearity at a given input value
If a nonlinear I/O system has no internal state, then evaluating the
system at an input `u` gives the output `y = F(u)`, determined by the
output function.
Parameters
----------
params : dict, optional
Parameter values for the system. Passed to the evaluation function
for the system as default values, overriding internal defaults.
squeeze : bool, optional
If True and if the system has a single output, return the system
output as a 1D array rather than a 2D array. If False, return the
system output as a 2D array even if the system is SISO. Default
value set by config.defaults['control.squeeze_time_response'].
"""
# Make sure the call makes sense
if not sys._isstatic():
raise TypeError(
"function evaluation is only supported for static "
"input/output systems")
# If we received any parameters, update them before calling _out()
if params is not None:
sys._update_params(params)
# Evaluate the function on the argument
out = sys._out(0, np.array((0,)), np.asarray(u))
_, out = _process_time_response(
None, out, issiso=sys.issiso(), squeeze=squeeze)
return out
def __mul__(self, other):
"""Multiply two input/output systems (series interconnection)"""
# Convert 'other' to an I/O system if needed
other = _convert_static_iosystem(other)
if not isinstance(other, InputOutputSystem):
return NotImplemented
# Make sure systems can be interconnected
if other.noutputs != self.ninputs:
raise ValueError(
"can't multiply systems with incompatible inputs and outputs")
# Make sure timebase are compatible
dt = common_timebase(other.dt, self.dt)
# Create a new system to handle the composition
inplist = [(0, i) for i in range(other.ninputs)]
outlist = [(1, i) for i in range(self.noutputs)]
newsys = InterconnectedSystem(
(other, self), inplist=inplist, outlist=outlist)
# Set up the connection map manually
newsys.set_connect_map(np.block(
[[np.zeros((other.ninputs, other.noutputs)),
np.zeros((other.ninputs, self.noutputs))],
[np.eye(self.ninputs, other.noutputs),
np.zeros((self.ninputs, self.noutputs))]]
))
# Return the newly created InterconnectedSystem
return newsys
def __rmul__(self, other):
"""Pre-multiply an input/output systems by a scalar/matrix"""
# Convert other to an I/O system if needed
other = _convert_static_iosystem(other)
if not isinstance(other, InputOutputSystem):
return NotImplemented
# Make sure systems can be interconnected
if self.noutputs != other.ninputs:
raise ValueError("Can't multiply systems with incompatible "
"inputs and outputs")
# Make sure timebase are compatible
dt = common_timebase(self.dt, other.dt)
# Create a new system to handle the composition
inplist = [(0, i) for i in range(self.ninputs)]
outlist = [(1, i) for i in range(other.noutputs)]
newsys = InterconnectedSystem(
(self, other), inplist=inplist, outlist=outlist)
# Set up the connection map manually
newsys.set_connect_map(np.block(
[[np.zeros((self.ninputs, self.noutputs)),
np.zeros((self.ninputs, other.noutputs))],
[np.eye(self.ninputs, self.noutputs),
np.zeros((other.ninputs, other.noutputs))]]
))
# Return the newly created InterconnectedSystem
return newsys
def __add__(self, other):
"""Add two input/output systems (parallel interconnection)"""
# Convert other to an I/O system if needed
other = _convert_static_iosystem(other)
if not isinstance(other, InputOutputSystem):
return NotImplemented
# Make sure number of input and outputs match
if self.ninputs != other.ninputs or self.noutputs != other.noutputs:
raise ValueError("Can't add systems with incompatible numbers of "
"inputs or outputs")
# Create a new system to handle the composition
inplist = [[(0, i), (1, i)] for i in range(self.ninputs)]
outlist = [[(0, i), (1, i)] for i in range(self.noutputs)]
newsys = InterconnectedSystem(
(self, other), inplist=inplist, outlist=outlist)
# Return the newly created InterconnectedSystem
return newsys
def __radd__(self, other):
"""Parallel addition of input/output system to a compatible object."""
# Convert other to an I/O system if needed
other = _convert_static_iosystem(other)
if not isinstance(other, InputOutputSystem):
return NotImplemented
# Make sure number of input and outputs match
if self.ninputs != other.ninputs or self.noutputs != other.noutputs:
raise ValueError("can't add systems with incompatible numbers of "
"inputs or outputs")
# Create a new system to handle the composition
inplist = [[(0, i), (1, i)] for i in range(other.ninputs)]
outlist = [[(0, i), (1, i)] for i in range(other.noutputs)]
newsys = InterconnectedSystem(
(other, self), inplist=inplist, outlist=outlist)
# Return the newly created InterconnectedSystem
return newsys
def __sub__(self, other):
"""Subtract two input/output systems (parallel interconnection)"""
# Convert other to an I/O system if needed
other = _convert_static_iosystem(other)
if not isinstance(other, InputOutputSystem):
return NotImplemented
# Make sure number of input and outputs match
if self.ninputs != other.ninputs or self.noutputs != other.noutputs:
raise ValueError(
"can't substract systems with incompatible numbers of "
"inputs or outputs")
ninputs = self.ninputs
noutputs = self.noutputs
# Create a new system to handle the composition
inplist = [[(0, i), (1, i)] for i in range(ninputs)]
outlist = [[(0, i), (1, i, -1)] for i in range(noutputs)]
newsys = InterconnectedSystem(
(self, other), inplist=inplist, outlist=outlist)
# Return the newly created InterconnectedSystem
return newsys
def __rsub__(self, other):
"""Parallel subtraction of I/O system to a compatible object."""
# Convert other to an I/O system if needed
other = _convert_static_iosystem(other)
if not isinstance(other, InputOutputSystem):
return NotImplemented
return other - self
def __neg__(self):
"""Negate an input/output system (rescale)"""
if self.ninputs is None or self.noutputs is None:
raise ValueError("Can't determine number of inputs or outputs")
# Create a new selftem to hold the negation
inplist = [(0, i) for i in range(self.ninputs)]
outlist = [(0, i, -1) for i in range(self.noutputs)]
newsys = InterconnectedSystem(
(self,), dt=self.dt, inplist=inplist, outlist=outlist)
# Return the newly created system
return newsys
def __truediv__(self, other):
"""Division of input/output system (by scalar or array)"""
if not isinstance(other, InputOutputSystem):
return self * (1/other)
else:
return NotImplemented
def _update_params(self, params, warning=False):
# Update the current parameter values
self._current_params = self.params.copy()
if params:
self._current_params.update(params)
def _rhs(self, t, x, u):
"""Evaluate right hand side of a differential or difference equation.
Private function used to compute the right hand side of an
input/output system model. Intended for fast evaluation; for a more
user-friendly interface you may want to use :meth:`dynamics`.
"""
return np.asarray(
self.updfcn(t, x, u, self._current_params)).reshape(-1)
def dynamics(self, t, x, u, params=None):
"""Compute the dynamics of a differential or difference equation.
Given time `t`, input `u` and state `x`, returns the value of the
right hand side of the dynamical system. If the system is continuous,
returns the time derivative
dx/dt = f(t, x, u[, params])
where `f` is the system's (possibly nonlinear) dynamics function.
If the system is discrete-time, returns the next value of `x`:
x[t+dt] = f(t, x[t], u[t][, params])
where `t` is a scalar.
The inputs `x` and `u` must be of the correct length. The `params`
argument is an optional dictionary of parameter values.
Parameters
----------
t : float
the time at which to evaluate
x : array_like
current state
u : array_like
input
params : dict, optional
system parameter values
Returns
-------
dx/dt or x[t+dt] : ndarray
"""
self._update_params(params)
return self._rhs(
t, np.asarray(x).reshape(-1), np.asarray(u).reshape(-1))
def _out(self, t, x, u):
"""Evaluate the output of a system at a given state, input, and time
Private function used to compute the output of of an input/output
system model given the state, input, parameters. Intended for fast
evaluation; for a more user-friendly interface you may want to use
:meth:`output`.
"""
#
# To allow lazy evaluation of the system size, we allow for the
# possibility that noutputs is left unspecified when the system
# is created => we have to check for that case here (and return
# the system state or a portion of it).
#
if self.outfcn is None:
return x if self.noutputs is None else x[:self.noutputs]
else:
return np.asarray(
self.outfcn(t, x, u, self._current_params)).reshape(-1)
def output(self, t, x, u, params=None):
"""Compute the output of the system
Given time `t`, input `u` and state `x`, returns the output of the
system:
y = g(t, x, u[, params])
The inputs `x` and `u` must be of the correct length.
Parameters
----------
t : float
the time at which to evaluate
x : array_like
current state
u : array_like
input
params : dict, optional
system parameter values
Returns
-------
y : ndarray
"""
self._update_params(params)
return self._out(
t, np.asarray(x).reshape(-1), np.asarray(u).reshape(-1))
def feedback(self, other=1, sign=-1, params=None):
"""Feedback interconnection between two input/output systems
Parameters
----------
sys1: InputOutputSystem
The primary process.
sys2: InputOutputSystem
The feedback process (often a feedback controller).
sign: scalar, optional
The sign of feedback. `sign` = -1 indicates negative feedback,
and `sign` = 1 indicates positive feedback. `sign` is an optional
argument; it assumes a value of -1 if not specified.
Returns
-------
out: InputOutputSystem
Raises
------
ValueError
if the inputs, outputs, or timebases of the systems are
incompatible.
"""
# Convert sys2 to an I/O system if needed
other = _convert_static_iosystem(other)
# Make sure systems can be interconnected
if self.noutputs != other.ninputs or other.noutputs != self.ninputs:
raise ValueError("Can't connect systems with incompatible "
"inputs and outputs")
# Make sure timebases are compatible
dt = common_timebase(self.dt, other.dt)
inplist = [(0, i) for i in range(self.ninputs)]
outlist = [(0, i) for i in range(self.noutputs)]
# Return the series interconnection between the systems
newsys = InterconnectedSystem(
(self, other), inplist=inplist, outlist=outlist,
params=params, dt=dt)
# Set up the connecton map manually
newsys.set_connect_map(np.block(
[[np.zeros((self.ninputs, self.noutputs)),
sign * np.eye(self.ninputs, other.noutputs)],
[np.eye(other.ninputs, self.noutputs),
np.zeros((other.ninputs, other.noutputs))]]
))
# Return the newly created system
return newsys
def linearize(self, x0, u0=None, t=0, params=None, eps=1e-6,
copy_names=False, **kwargs):
"""Linearize an input/output system at a given state and input.
Return the linearization of an input/output system at a given
operating point (or state and input value) as a StateSpace system.
See :func:`~control.linearize` for complete documentation.
"""
#
# Default method: if the linearization is not defined by the
# subclass, perform a numerical linearization use the `_rhs()` and
# `_out()` member functions.
#
from .statesp import StateSpace
# Allow first argument to be an operating point
if isinstance(x0, OperatingPoint):
u0 = x0.inputs if u0 is None else u0
x0 = x0.states
elif u0 is None:
u0 = 0
# Process nominal states and inputs
x0, nstates = _process_vector_argument(x0, "x0", self.nstates)
u0, ninputs = _process_vector_argument(u0, "u0", self.ninputs)
# Update the current parameters (prior to calling _out())
self._update_params(params)
# Compute number of outputs by evaluating the output function
noutputs = _find_size(self.noutputs, self._out(t, x0, u0), "outputs")
# Compute the nominal value of the update law and output
F0 = self._rhs(t, x0, u0)
H0 = self._out(t, x0, u0)
# Create empty matrices that we can fill up with linearizations
A = np.zeros((nstates, nstates)) # Dynamics matrix
B = np.zeros((nstates, ninputs)) # Input matrix
C = np.zeros((noutputs, nstates)) # Output matrix
D = np.zeros((noutputs, ninputs)) # Direct term
# Perturb each of the state variables and compute linearization
for i in range(nstates):
dx = np.zeros((nstates,))
dx[i] = eps
A[:, i] = (self._rhs(t, x0 + dx, u0) - F0) / eps
C[:, i] = (self._out(t, x0 + dx, u0) - H0) / eps
# Perturb each of the input variables and compute linearization
for i in range(ninputs):
du = np.zeros((ninputs,))
du[i] = eps
B[:, i] = (self._rhs(t, x0, u0 + du) - F0) / eps
D[:, i] = (self._out(t, x0, u0 + du) - H0) / eps
# Create the state space system
linsys = StateSpace(A, B, C, D, self.dt, remove_useless_states=False)
# Set the system name, inputs, outputs, and states
if copy_names:
linsys._copy_names(self, prefix_suffix_name='linearized')
# re-init to include desired signal names if names were provided
return StateSpace(linsys, **kwargs)
class InterconnectedSystem(NonlinearIOSystem):
"""Interconnection of a set of input/output systems.
This class is used to implement a system that is an interconnection of
input/output systems. The sys consists of a collection of subsystems
whose inputs and outputs are connected via a connection map. The overall
system inputs and outputs are subsets of the subsystem inputs and outputs.
The function :func:`~control.interconnect` should be used to create an
interconnected I/O system since it performs additional argument
processing and checking.
"""
def __init__(self, syslist, connections=None, inplist=None, outlist=None,
params=None, warn_duplicate=None, connection_type=None,
**kwargs):
"""Create an I/O system from a list of systems + connection info."""
from .statesp import _convert_to_statespace
from .xferfcn import TransferFunction
self.connection_type = connection_type # explicit, implicit, or None
# Convert input and output names to lists if they aren't already
if inplist is not None and not isinstance(inplist, list):
inplist = [inplist]
if outlist is not None and not isinstance(outlist, list):
outlist = [outlist]
# Check if dt argument was given; if not, pull from systems
dt = kwargs.pop('dt', None)
# Process keyword arguments (except dt)
name, inputs, outputs, states, _ = _process_iosys_keywords(kwargs)
# Initialize the system list and index
self.syslist = list(syslist) # ensure modifications can be made
self.syslist_index = {}
# Initialize the input, output, and state counts, indices
nstates, self.state_offset = 0, []
ninputs, self.input_offset = 0, []
noutputs, self.output_offset = 0, []
# Keep track of system objects and names we have already seen
sysobj_name_dct = {}
sysname_count_dct = {}
# Go through the system list and keep track of counts, offsets
for sysidx, sys in enumerate(self.syslist):
# Convert transfer functions to state space
if isinstance(sys, TransferFunction):
sys = _convert_to_statespace(sys)
self.syslist[sysidx] = sys
# Make sure time bases are consistent
dt = common_timebase(dt, sys.dt)
# Make sure number of inputs, outputs, states is given
if sys.ninputs is None or sys.noutputs is None:
raise TypeError("system '%s' must define number of inputs, "
"outputs, states in order to be connected" %
sys.name)
elif sys.nstates is None:
raise TypeError("can't interconnect systems with no state")
# Keep track of the offsets into the states, inputs, outputs
self.input_offset.append(ninputs)
self.output_offset.append(noutputs)
self.state_offset.append(nstates)
# Keep track of the total number of states, inputs, outputs
nstates += sys.nstates
ninputs += sys.ninputs
noutputs += sys.noutputs
# Check for duplicate systems or duplicate names
# Duplicates are renamed sysname_1, sysname_2, etc.
if sys in sysobj_name_dct:
# Make a copy of the object using a new name
if warn_duplicate is None and sys._generic_name_check():
# Make a copy w/out warning, using generic format
sys = sys.copy(use_prefix_suffix=False)
warn_flag = False
else:
sys = sys.copy()
warn_flag = warn_duplicate
# Warn the user about the new object
if warn_flag is not False:
warn("duplicate object found in system list; "
"created copy: %s" % str(sys.name), stacklevel=2)
# Check to see if the system name shows up more than once
if sys.name is not None and sys.name in sysname_count_dct:
count = sysname_count_dct[sys.name]
sysname_count_dct[sys.name] += 1
sysname = sys.name + "_" + str(count)
sysobj_name_dct[sys] = sysname
self.syslist_index[sysname] = sysidx
if warn_duplicate is not False:
warn("duplicate name found in system list; "
"renamed to {}".format(sysname), stacklevel=2)
else:
sysname_count_dct[sys.name] = 1
sysobj_name_dct[sys] = sys.name
self.syslist_index[sys.name] = sysidx
if states is None:
states = []
state_name_delim = config.defaults['iosys.state_name_delim']
for sys, sysname in sysobj_name_dct.items():
states += [sysname + state_name_delim +
statename for statename in sys.state_index.keys()]
# Make sure we the state list is the right length (internal check)
if isinstance(states, list) and len(states) != nstates:
raise RuntimeError(
f"construction of state labels failed; found: "
f"{len(states)} labels; expecting {nstates}")
# Figure out what the inputs and outputs are
if inputs is None and inplist is not None:
inputs = len(inplist)
if outputs is None and outlist is not None:
outputs = len(outlist)
# Create updfcn and outfcn
def updfcn(t, x, u, params):
self._update_params(params)
return self._rhs(t, x, u)
def outfcn(t, x, u, params):
self._update_params(params)
return self._out(t, x, u)
# Initialize NonlinearIOSystem object
super().__init__(
updfcn, outfcn, inputs=inputs, outputs=outputs,
states=states, dt=dt, name=name, params=params, **kwargs)
# Convert the list of interconnections to a connection map (matrix)
self.connect_map = np.zeros((ninputs, noutputs))
for connection in connections or []:
input_indices = self._parse_input_spec(connection[0])
for output_spec in connection[1:]:
output_indices, gain = self._parse_output_spec(output_spec)
if len(output_indices) != len(input_indices):
raise ValueError(
f"inconsistent number of signals in connecting"
f" '{output_spec}' to '{connection[0]}'")
for input_index, output_index in zip(
input_indices, output_indices):
if self.connect_map[input_index, output_index] != 0:
warn("multiple connections given for input %d" %
input_index + "; combining with previous entries")
self.connect_map[input_index, output_index] += gain
# Convert the input list to a matrix: maps system to subsystems
self.input_map = np.zeros((ninputs, self.ninputs))
for index, inpspec in enumerate(inplist or []):
if isinstance(inpspec, (int, str, tuple)):
inpspec = [inpspec]
if not isinstance(inpspec, list):
raise ValueError("specifications in inplist must be of type "
"int, str, tuple or list")
for spec in inpspec:
ulist_indices = self._parse_input_spec(spec)
for j, ulist_index in enumerate(ulist_indices):
if self.input_map[ulist_index, index] != 0:
warn("multiple connections given for input %d" %
index + "; combining with previous entries.")
self.input_map[ulist_index, index + j] += 1
# Convert the output list to a matrix: maps subsystems to system
self.output_map = np.zeros((self.noutputs, noutputs + ninputs))
for index, outspec in enumerate(outlist or []):
if isinstance(outspec, (int, str, tuple)):
outspec = [outspec]
if not isinstance(outspec, list):
raise ValueError("specifications in outlist must be of type "
"int, str, tuple or list")
for spec in outspec:
ylist_indices, gain = self._parse_output_spec(spec)
for j, ylist_index in enumerate(ylist_indices):
if self.output_map[index, ylist_index] != 0:
warn("multiple connections given for output %d" %
index + "; combining with previous entries")
self.output_map[index + j, ylist_index] += gain
def _update_params(self, params, warning=False):
for sys in self.syslist:
local = sys.params.copy() # start with system parameters
local.update(self.params) # update with global params
if params:
local.update(params) # update with locally passed parameters
sys._update_params(local, warning=warning)
def _rhs(self, t, x, u):
# Make sure state and input are vectors
x = np.array(x, ndmin=1)
u = np.array(u, ndmin=1)
# Compute the input and output vectors
ulist, ylist = self._compute_static_io(t, x, u)
# Go through each system and update the right hand side for that system
xdot = np.zeros((self.nstates,)) # Array to hold results
state_index, input_index = 0, 0 # Start at the beginning
for sys in self.syslist:
# Update the right hand side for this subsystem
if sys.nstates != 0:
xdot[state_index:state_index + sys.nstates] = sys._rhs(
t, x[state_index:state_index + sys.nstates],
ulist[input_index:input_index + sys.ninputs])
# Update the state and input index counters
state_index += sys.nstates
input_index += sys.ninputs
return xdot
def _out(self, t, x, u):
# Make sure state and input are vectors
x = np.array(x, ndmin=1)
u = np.array(u, ndmin=1)
# Compute the input and output vectors
ulist, ylist = self._compute_static_io(t, x, u)
# Make the full set of subsystem outputs to system output
return self.output_map @ ylist
def _compute_static_io(self, t, x, u):
# Figure out the total number of inputs and outputs
(ninputs, noutputs) = self.connect_map.shape
#
# Get the outputs and inputs at the current system state
#
# Initialize the lists used to keep track of internal signals
ulist = np.dot(self.input_map, u)
ylist = np.zeros((noutputs + ninputs,))
# To allow for feedthrough terms, iterate multiple times to allow
# feedthrough elements to propagate. For n systems, we could need to
# cycle through n+1 times before reaching steady state
# TODO (later): see if there is a more efficient way to compute
cycle_count = len(self.syslist) + 1
while cycle_count > 0:
state_index, input_index, output_index = 0, 0, 0
for sys in self.syslist:
# Compute outputs for each system from current state
ysys = sys._out(
t, x[state_index:state_index + sys.nstates],
ulist[input_index:input_index + sys.ninputs])
# Store the outputs at the start of ylist
ylist[output_index:output_index + sys.noutputs] = \
ysys.reshape((-1,))
# Store the input in the second part of ylist
ylist[noutputs + input_index:
noutputs + input_index + sys.ninputs] = \
ulist[input_index:input_index + sys.ninputs]
# Increment the index pointers
state_index += sys.nstates
input_index += sys.ninputs
output_index += sys.noutputs
# Compute inputs based on connection map
new_ulist = self.connect_map @ ylist[:noutputs] \
+ np.dot(self.input_map, u)
# Check to see if any of the inputs changed
if (ulist == new_ulist).all():
break
else:
ulist = new_ulist
# Decrease the cycle counter
cycle_count -= 1
# Make sure that we stopped before detecting an algebraic loop
if cycle_count == 0:
raise RuntimeError("algebraic loop detected")
return ulist, ylist
def _parse_input_spec(self, spec):
"""Parse an input specification and returns the indices."""
# Parse the signal that we received
subsys_index, input_indices, gain = _parse_spec(
self.syslist, spec, 'input')
if gain != 1:
raise ValueError("gain not allowed in spec '%s'" % str(spec))
# Return the indices into the input vector list (ylist)
return [self.input_offset[subsys_index] + i for i in input_indices]
def _parse_output_spec(self, spec):
"""Parse an output specification and returns the indices and gain."""
# Parse the rest of the spec with standard signal parsing routine
try:
# Start by looking in the set of subsystem outputs
subsys_index, output_indices, gain = \
_parse_spec(self.syslist, spec, 'output')
output_offset = self.output_offset[subsys_index]
except ValueError:
# Try looking in the set of subsystem *inputs*
subsys_index, output_indices, gain = _parse_spec(
self.syslist, spec, 'input or output', dictname='input_index')
# Return the index into the input vector list (ylist)
output_offset = sum(sys.noutputs for sys in self.syslist) + \
self.input_offset[subsys_index]
return [output_offset + i for i in output_indices], gain
def _find_system(self, name):
return self.syslist_index.get(name, None)
def set_connect_map(self, connect_map):
"""Set the connection map for an interconnected I/O system.
Parameters
----------
connect_map : 2D array
Specify the matrix that will be used to multiply the vector of
subsystem outputs to obtain the vector of subsystem inputs.
"""
# Make sure the connection map is the right size
if connect_map.shape != self.connect_map.shape:
ValueError("Connection map is not the right shape")
self.connect_map = connect_map
def set_input_map(self, input_map):
"""Set the input map for an interconnected I/O system.
Parameters
----------
input_map : 2D array
Specify the matrix that will be used to multiply the vector of
system inputs to obtain the vector of subsystem inputs. These
values are added to the inputs specified in the connection map.
"""
# Figure out the number of internal inputs
ninputs = sum(sys.ninputs for sys in self.syslist)
# Make sure the input map is the right size
if input_map.shape[0] != ninputs:
ValueError("Input map is not the right shape")
self.input_map = input_map
self.ninputs = input_map.shape[1]
def set_output_map(self, output_map):
"""Set the output map for an interconnected I/O system.
Parameters
----------
output_map : 2D array
Specify the matrix that will be used to multiply the vector of
subsystem outputs concatenated with subsystem inputs to obtain
the vector of system outputs.
"""
# Figure out the number of internal inputs and outputs
ninputs = sum(sys.ninputs for sys in self.syslist)
noutputs = sum(sys.noutputs for sys in self.syslist)
# Make sure the output map is the right size
if output_map.shape[1] == noutputs:
# For backward compatibility, add zeros to the end of the array
output_map = np.concatenate(
(output_map,
np.zeros((output_map.shape[0], ninputs))),
axis=1)
if output_map.shape[1] != noutputs + ninputs:
ValueError("Output map is not the right shape")
self.output_map = output_map
self.noutputs = output_map.shape[0]
def unused_signals(self):
"""Find unused subsystem inputs and outputs
Returns
-------
unused_inputs : dict
A mapping from tuple of indices (isys, isig) to string
'{sys}.{sig}', for all unused subsystem inputs.
unused_outputs : dict
A mapping from tuple of indices (osys, osig) to string
'{sys}.{sig}', for all unused subsystem outputs.
"""
used_sysinp_via_inp = np.nonzero(self.input_map)[0]
used_sysout_via_out = np.nonzero(self.output_map)[1]
used_sysinp_via_con, used_sysout_via_con = np.nonzero(self.connect_map)
used_sysinp = set(used_sysinp_via_inp) | set(used_sysinp_via_con)
used_sysout = set(used_sysout_via_out) | set(used_sysout_via_con)
nsubsysinp = sum(sys.ninputs for sys in self.syslist)
nsubsysout = sum(sys.noutputs for sys in self.syslist)