python-control/control/xferfcn.py at f98f9d76d9f051dd96539260ad3658bdfa25190f · python-control/python-control

Repository navigation

History

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

"""xferfcn.py

Transfer function representation and functions.

This file contains the TransferFunction class and also functions

that operate on transfer functions. This is the primary representation

for the python-control library.

"""

Redistribution and use in source and binary forms, with or without

modification, are permitted provided that the following conditions

are met:

1. Redistributions of source code must retain the above copyright

notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright

notice, this list of conditions and the following disclaimer in the

documentation and/or other materials provided with the distribution.

3. Neither the name of the California Institute of Technology nor

the names of its contributors may be used to endorse or promote

products derived from this software without specific prior

written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS

"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT

LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS

FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CALTECH

OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,

SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT

LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF

USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND

ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,

OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT

OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF

SUCH DAMAGE.

Author: Richard M. Murray

Date: 24 May 09

Revised: Kevin K. Chen, Dec 10

$Id$

"""

# External function declarations

import numpy as np

from numpy import angle, array, empty, finfo, ndarray, ones, \

polyadd, polymul, polyval, roots, sqrt, zeros, squeeze, exp, pi, \

where, delete, real, poly, nonzero

import scipy as sp

from scipy.signal import tf2zpk, zpk2tf, cont2discrete

from scipy.signal import TransferFunction as signalTransferFunction

from copy import deepcopy

from warnings import warn

from itertools import chain

from re import sub

from .lti import LTI, _process_frequency_response

from .iosys import InputOutputSystem, common_timebase, isdtime, \

_process_iosys_keywords

from .exception import ControlMIMONotImplemented

from .frdata import FrequencyResponseData

from . import config

__all__ = ['TransferFunction', 'tf', 'zpk', 'ss2tf', 'tfdata']

# Define module default parameter values

_xferfcn_defaults = {

'xferfcn.display_format': 'poly',

'xferfcn.floating_point_format': '.4g'

}

class TransferFunction(LTI):

"""TransferFunction(num, den[, dt])

A class for representing transfer functions.

The TransferFunction class is used to represent systems in transfer

function form.

Parameters

----------

num : array_like, or list of list of array_like

Polynomial coefficients of the numerator

den : array_like, or list of list of array_like

Polynomial coefficients of the denominator

dt : None, True or float, optional

System timebase. 0 (default) indicates continuous

time, True indicates discrete time with unspecified sampling

time, positive number is discrete time with specified

sampling time, None indicates unspecified timebase (either

continuous or discrete time).

display_format: None, 'poly' or 'zpk'

Set the display format used in printing the TransferFunction object.

Default behavior is polynomial display and can be changed by

changing config.defaults['xferfcn.display_format'].

Attributes

----------

ninputs, noutputs, nstates : int

Number of input, output and state variables.

num, den : 2D list of array

Polynomial coefficients of the numerator and denominator.

dt : None, True or float

System timebase. 0 (default) indicates continuous time, True indicates

discrete time with unspecified sampling time, positive number is

discrete time with specified sampling time, None indicates unspecified

timebase (either continuous or discrete time).

Notes

-----

The attribues 'num' and 'den' are 2-D lists of arrays containing MIMO

numerator and denominator coefficients. For example,

>>> num[2][5] = numpy.array([1., 4., 8.]) # doctest: +SKIP

means that the numerator of the transfer function from the 6th input to

the 3rd output is set to s^2 + 4s + 8.

A discrete time transfer function is created by specifying a nonzero

'timebase' dt when the system is constructed:

* 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

Systems must have compatible timebases in order to be combined. A discrete

time system with unspecified sampling time (`dt = True`) can be combined

with a system having a specified sampling time; the result will be a

discrete time system with the sample time of the latter system. Similarly,

a system with timebase `None` can be combined with a system having any

timebase; the result will have the timebase of the latter system.

The default value of dt can be changed by changing the value of

``control.config.defaults['control.default_dt']``.

A transfer function is callable and returns the value of the transfer

function evaluated at a point in the complex plane. See

:meth:`~control.TransferFunction.__call__` for a more detailed description.

The TransferFunction class defines two constants ``s`` and ``z`` that

represent the differentiation and delay operators in continuous and

discrete time. These can be used to create variables that allow algebraic

creation of transfer functions. For example,

>>> s = ct.TransferFunction.s

>>> G = (s + 1)/(s**2 + 2*s + 1)

"""

def __init__(self, *args, **kwargs):

"""TransferFunction(num, den[, dt])

Construct a transfer function.

The default constructor is TransferFunction(num, den), where num and

den are lists of lists of arrays containing polynomial coefficients.

To create a discrete time transfer funtion, use TransferFunction(num,

den, dt) where 'dt' is the sampling time (or True for unspecified

sampling time). To call the copy constructor, call

TransferFunction(sys), where sys is a TransferFunction object

(continuous or discrete).

"""

# Process positional arguments

if len(args) == 2:

# The user provided a numerator and a denominator.

num, den = args

elif len(args) == 3:

# Discrete time transfer function

num, den, dt = args

if 'dt' in kwargs:

warn("received multiple dt arguments, "

"using positional arg dt = %s" % dt)

kwargs['dt'] = dt

args = args[:-1]

elif len(args) == 1:

# Use the copy constructor.

if not isinstance(args[0], TransferFunction):

raise TypeError("The one-argument constructor can only take \

in a TransferFunction object. Received %s."

% type(args[0]))

num = args[0].num

den = args[0].den

else:

raise TypeError("Needs 1, 2 or 3 arguments; received %i."

% len(args))

num = _clean_part(num)

den = _clean_part(den)

# Process keyword arguments

# During module init, TransferFunction.s and TransferFunction.z

# get initialized when defaults are not fully initialized yet.

# Use 'poly' in these cases.

self.display_format = kwargs.pop(

'display_format',

config.defaults.get('xferfcn.display_format', 'poly'))

if self.display_format not in ('poly', 'zpk'):

raise ValueError("display_format must be 'poly' or 'zpk',"

" got '%s'" % self.display_format)

# Determine if the transfer function is static (needed for dt)

static = True

for col in num + den:

for poly in col:

if len(poly) > 1:

static = False

defaults = args[0] if len(args) == 1 else \

{'inputs': len(num[0]), 'outputs': len(num)}

name, inputs, outputs, states, dt = _process_iosys_keywords(

kwargs, defaults, static=static)

if states:

raise TypeError(

"states keyword not allowed for transfer functions")

# Initialize LTI (InputOutputSystem) object

super().__init__(

name=name, inputs=inputs, outputs=outputs, dt=dt, **kwargs)

# Check to make sure everything is consistent

# Make sure numerator and denominator matrices have consistent sizes

if self.ninputs != len(den[0]):

raise ValueError(

"The numerator has %i input(s), but the denominator has "

"%i input(s)." % (self.ninputs, len(den[0])))

if self.noutputs != len(den):

raise ValueError(

"The numerator has %i output(s), but the denominator has "

"%i output(s)." % (self.noutputs, len(den)))

# Additional checks/updates on structure of the transfer function

for i in range(self.noutputs):

# Make sure that each row has the same number of columns

if len(num[i]) != self.ninputs:

raise ValueError(

"Row 0 of the numerator matrix has %i elements, but row "

"%i has %i." % (self.ninputs, i, len(num[i])))

if len(den[i]) != self.ninputs:

raise ValueError(

"Row 0 of the denominator matrix has %i elements, but row "

"%i has %i." % (self.ninputs, i, len(den[i])))

# Check for zeros in numerator or denominator

# TODO: Right now these checks are only done during construction.

# It might be worthwhile to think of a way to perform checks if the

# user modifies the transfer function after construction.

for j in range(self.ninputs):

# Check that we don't have any zero denominators.

zeroden = True

for k in den[i][j]:

if k:

zeroden = False

break

if zeroden:

raise ValueError(

"Input %i, output %i has a zero denominator."

% (j + 1, i + 1))

# If we have zero numerators, set the denominator to 1.

zeronum = True

for k in num[i][j]:

if k:

zeronum = False

break

if zeronum:

den[i][j] = ones(1)

# Store the numerator and denominator

self.num = num

self.den = den

# Final processing

# Truncate leading zeros

self._truncatecoeff()

# Class attributes

# These attributes are defined as class attributes so that they are

# documented properly. They are "overwritten" in __init__.

#: Number of system inputs.

#: :meta hide-value:

ninputs = 1

#: Number of system outputs.

#: :meta hide-value:

noutputs = 1

#: Transfer function numerator polynomial (array)

#: The numerator of the transfer function is stored as an 2D list of

#: arrays containing MIMO numerator coefficients, indexed by outputs and

#: inputs. For example, ``num[2][5]`` is the array of coefficients for

#: the numerator of the transfer function from the sixth input to the

#: third output.

#: :meta hide-value:

num = [[0]]

#: Transfer function denominator polynomial (array)

#: The numerator of the transfer function is store as an 2D list of

#: arrays containing MIMO numerator coefficients, indexed by outputs and

#: inputs. For example, ``den[2][5]`` is the array of coefficients for

#: the denominator of the transfer function from the sixth input to the

#: third output.

#: :meta hide-value:

den = [[0]]

def __call__(self, x, squeeze=None, warn_infinite=True):

"""Evaluate system's transfer function at complex frequencies.

Returns the complex frequency response `sys(x)` where `x` is `s` for

continuous-time systems and `z` for discrete-time systems.

In general the system may be multiple input, multiple output

(MIMO), where `m = self.ninputs` number of inputs and `p =

self.noutputs` number of outputs.

To evaluate at a frequency omega in radians per second, enter

``x = omega * 1j``, for continuous-time systems, or

``x = exp(1j * omega * dt)`` for discrete-time systems. Or use

:meth:`TransferFunction.frequency_response`.

Parameters

----------

x : complex or complex 1D array_like

Complex frequencies

squeeze : bool, optional

If squeeze=True, remove single-dimensional entries from the shape

of the output even if the system is not SISO. If squeeze=False,

keep all indices (output, input and, if omega is array_like,

frequency) even if the system is SISO. The default value can be

set using config.defaults['control.squeeze_frequency_response'].

If True and the system is single-input single-output (SISO),

return a 1D array rather than a 3D array. Default value (True)

set by config.defaults['control.squeeze_frequency_response'].

warn_infinite : bool, optional

If set to `False`, turn off divide by zero warning.

Returns

-------

fresp : complex ndarray

The frequency response of the system. If the system is SISO and

squeeze is not True, the shape of the array matches the shape of

omega. If the system is not SISO or squeeze is False, the first

two dimensions of the array are indices for the output and input

and the remaining dimensions match omega. If ``squeeze`` is True

then single-dimensional axes are removed.

"""

out = self.horner(x, warn_infinite=warn_infinite)

return _process_frequency_response(self, x, out, squeeze=squeeze)

def horner(self, x, warn_infinite=True):

"""Evaluate system's transfer function at complex frequency

using Horner's method.

Evaluates `sys(x)` where `x` is `s` for continuous-time systems and `z`

for discrete-time systems.

Expects inputs and outputs to be formatted correctly. Use ``sys(x)``

for a more user-friendly interface.

Parameters

----------

x : complex array_like or complex scalar

Complex frequencies

Returns

-------

output : (self.noutputs, self.ninputs, len(x)) complex ndarray

Frequency response

"""

# Make sure the argument is a 1D array of complex numbers

x_arr = np.atleast_1d(x).astype(complex, copy=False)

# Make sure that we are operating on a simple list

if len(x_arr.shape) > 1:

raise ValueError("input list must be 1D")

# Initialize the output matrix in the proper shape

out = empty((self.noutputs, self.ninputs, len(x_arr)), dtype=complex)

# Set up error processing based on warn_infinite flag

with np.errstate(all='warn' if warn_infinite else 'ignore'):

for i in range(self.noutputs):

for j in range(self.ninputs):

out[i][j] = (polyval(self.num[i][j], x_arr) /

polyval(self.den[i][j], x_arr))

return out

def _truncatecoeff(self):

"""Remove extraneous zero coefficients from num and den.

Check every element of the numerator and denominator matrices, and

truncate leading zeros. For instance, running self._truncatecoeff()

will reduce self.num = [[[0, 0, 1, 2]]] to [[[1, 2]]].

"""

# Beware: this is a shallow copy. This should be okay.

data = [self.num, self.den]

for p in range(len(data)):

for i in range(self.noutputs):

for j in range(self.ninputs):

# Find the first nontrivial coefficient.

nonzero = None

for k in range(data[p][i][j].size):

if data[p][i][j][k]:

nonzero = k

break

if nonzero is None:

# The array is all zeros.

data[p][i][j] = zeros(1)

else:

# Truncate the trivial coefficients.

data[p][i][j] = data[p][i][j][nonzero:]

[self.num, self.den] = data

def __str__(self, var=None):

"""String representation of the transfer function.

Based on the display_format property, the output will be formatted as

either polynomials or in zpk form.

"""

mimo = not self.issiso()

if var is None:

var = 's' if self.isctime() else 'z'

outstr = f"{InputOutputSystem.__str__(self)}\n"

for ni in range(self.ninputs):

for no in range(self.noutputs):

if mimo:

outstr += "\nInput %i to output %i:" % (ni + 1, no + 1)

# Convert the numerator and denominator polynomials to strings.

if self.display_format == 'poly':

numstr = _tf_polynomial_to_string(self.num[no][ni], var=var)

denstr = _tf_polynomial_to_string(self.den[no][ni], var=var)

elif self.display_format == 'zpk':

num = self.num[no][ni]

if num.size == 1 and num.item() == 0:

# Catch a special case that SciPy doesn't handle

z, p, k = tf2zpk([1.], self.den[no][ni])

k = 0

else:

z, p, k = tf2zpk(self.num[no][ni], self.den[no][ni])

numstr = _tf_factorized_polynomial_to_string(

z, gain=k, var=var)

denstr = _tf_factorized_polynomial_to_string(p, var=var)

# Figure out the length of the separating line

dashcount = max(len(numstr), len(denstr))

dashes = '-' * dashcount

# Center the numerator or denominator

if len(numstr) < dashcount:

numstr = ' ' * ((dashcount - len(numstr)) // 2) + numstr

if len(denstr) < dashcount:

denstr = ' ' * ((dashcount - len(denstr)) // 2) + denstr

outstr += "\n" + numstr + "\n" + dashes + "\n" + denstr + "\n"

# If this is a strict discrete time system, print the sampling time

if type(self.dt) != bool and self.isdtime(strict=True):

outstr += "\ndt = " + str(self.dt) + "\n"

return outstr

# represent to implement a re-loadable version

def __repr__(self):

"""Print transfer function in loadable form"""

if self.issiso():

return "TransferFunction({num}, {den}{dt})".format(

num=self.num[0][0].__repr__(), den=self.den[0][0].__repr__(),

dt=', {}'.format(self.dt) if isdtime(self, strict=True)

else '')

else:

return "TransferFunction({num}, {den}{dt})".format(

num=self.num.__repr__(), den=self.den.__repr__(),

dt=', {}'.format(self.dt) if isdtime(self, strict=True)

else '')

def _repr_latex_(self, var=None):

"""LaTeX representation of transfer function, for Jupyter notebook"""

mimo = not self.issiso()

if var is None:

var = 's' if self.isctime() else 'z'

out = ['$$']

if mimo:

out.append(r"\begin{bmatrix}")

for no in range(self.noutputs):

for ni in range(self.ninputs):

# Convert the numerator and denominator polynomials to strings.

if self.display_format == 'poly':

numstr = _tf_polynomial_to_string(self.num[no][ni], var=var)

denstr = _tf_polynomial_to_string(self.den[no][ni], var=var)

elif self.display_format == 'zpk':

z, p, k = tf2zpk(self.num[no][ni], self.den[no][ni])

numstr = _tf_factorized_polynomial_to_string(

z, gain=k, var=var)

denstr = _tf_factorized_polynomial_to_string(p, var=var)

numstr = _tf_string_to_latex(numstr, var=var)

denstr = _tf_string_to_latex(denstr, var=var)

out += [r"\frac{", numstr, "}{", denstr, "}"]

if mimo and ni < self.ninputs - 1:

out.append("&")

if mimo:

out.append(r"\\")

if mimo:

out.append(r" \end{bmatrix}")

# See if this is a discrete time system with specific sampling time

if not (self.dt is None) and type(self.dt) != bool and self.dt > 0:

out += [r"\quad dt = ", str(self.dt)]

out.append("$$")

return ''.join(out)

def __neg__(self):

"""Negate a transfer function."""

num = deepcopy(self.num)

for i in range(self.noutputs):

for j in range(self.ninputs):

num[i][j] *= -1

return TransferFunction(num, self.den, self.dt)

def __add__(self, other):

"""Add two LTI objects (parallel connection)."""

from .statesp import StateSpace

# Convert the second argument to a transfer function.

if isinstance(other, StateSpace):

other = _convert_to_transfer_function(other)

elif isinstance(other, (int, float, complex, np.number, np.ndarray)):

other = _convert_to_transfer_function(other, inputs=self.ninputs,

outputs=self.noutputs)

if not isinstance(other, TransferFunction):

return NotImplemented

# Check that the input-output sizes are consistent.

if self.ninputs != other.ninputs:

raise ValueError(

"The first summand has %i input(s), but the second has %i."

% (self.ninputs, other.ninputs))

if self.noutputs != other.noutputs:

raise ValueError(

"The first summand has %i output(s), but the second has %i."

% (self.noutputs, other.noutputs))

dt = common_timebase(self.dt, other.dt)

# Preallocate the numerator and denominator of the sum.

num = [[[] for j in range(self.ninputs)] for i in range(self.noutputs)]

den = [[[] for j in range(self.ninputs)] for i in range(self.noutputs)]

for i in range(self.noutputs):

for j in range(self.ninputs):

num[i][j], den[i][j] = _add_siso(

self.num[i][j], self.den[i][j],

other.num[i][j], other.den[i][j])

return TransferFunction(num, den, dt)

def __radd__(self, other):

"""Right add two LTI objects (parallel connection)."""

return self + other

def __sub__(self, other):

"""Subtract two LTI objects."""

return self + (-other)

def __rsub__(self, other):

"""Right subtract two LTI objects."""

return other + (-self)

def __mul__(self, other):

"""Multiply two LTI objects (serial connection)."""

from .statesp import StateSpace

# Convert the second argument to a transfer function.

if isinstance(other, StateSpace):

other = _convert_to_transfer_function(other)

elif isinstance(other, (int, float, complex, np.number, np.ndarray)):

other = _convert_to_transfer_function(other, inputs=self.ninputs,

outputs=self.noutputs)

if not isinstance(other, TransferFunction):

return NotImplemented

# Check that the input-output sizes are consistent.

if self.ninputs != other.noutputs:

raise ValueError(

"C = A * B: A has %i column(s) (input(s)), but B has %i "

"row(s)\n(output(s))." % (self.ninputs, other.noutputs))

inputs = other.ninputs

outputs = self.noutputs

dt = common_timebase(self.dt, other.dt)

# Preallocate the numerator and denominator of the sum.

num = [[[0] for j in range(inputs)] for i in range(outputs)]

den = [[[1] for j in range(inputs)] for i in range(outputs)]

# Temporary storage for the summands needed to find the (i, j)th

# element of the product.

num_summand = [[] for k in range(self.ninputs)]

den_summand = [[] for k in range(self.ninputs)]

# Multiply & add.

for row in range(outputs):

for col in range(inputs):

for k in range(self.ninputs):

num_summand[k] = polymul(

self.num[row][k], other.num[k][col])

den_summand[k] = polymul(

self.den[row][k], other.den[k][col])

num[row][col], den[row][col] = _add_siso(

num[row][col], den[row][col],

num_summand[k], den_summand[k])

return TransferFunction(num, den, dt)

def __rmul__(self, other):

"""Right multiply two LTI objects (serial connection)."""

# Convert the second argument to a transfer function.

if isinstance(other, (int, float, complex, np.number)):

other = _convert_to_transfer_function(other, inputs=self.ninputs,

outputs=self.ninputs)

else:

other = _convert_to_transfer_function(other)

# Check that the input-output sizes are consistent.

if other.ninputs != self.noutputs:

raise ValueError(

"C = A * B: A has %i column(s) (input(s)), but B has %i "

"row(s)\n(output(s))." % (other.ninputs, self.noutputs))

inputs = self.ninputs

outputs = other.noutputs

dt = common_timebase(self.dt, other.dt)

# Preallocate the numerator and denominator of the sum.

num = [[[0] for j in range(inputs)] for i in range(outputs)]

den = [[[1] for j in range(inputs)] for i in range(outputs)]

# Temporary storage for the summands needed to find the

# (i, j)th element

# of the product.

num_summand = [[] for k in range(other.ninputs)]

den_summand = [[] for k in range(other.ninputs)]

for i in range(outputs): # Iterate through rows of product.

for j in range(inputs): # Iterate through columns of product.

for k in range(other.ninputs): # Multiply & add.

num_summand[k] = polymul(other.num[i][k], self.num[k][j])

den_summand[k] = polymul(other.den[i][k], self.den[k][j])

num[i][j], den[i][j] = _add_siso(

num[i][j], den[i][j],

num_summand[k], den_summand[k])

return TransferFunction(num, den, dt)

# TODO: Division of MIMO transfer function objects is not written yet.

def __truediv__(self, other):

"""Divide two LTI objects."""

if isinstance(other, (int, float, complex, np.number)):

other = _convert_to_transfer_function(

other, inputs=self.ninputs,

outputs=self.ninputs)

else:

other = _convert_to_transfer_function(other)

if (self.ninputs > 1 or self.noutputs > 1 or

other.ninputs > 1 or other.noutputs > 1):

raise NotImplementedError(

"TransferFunction.__truediv__ is currently \

implemented only for SISO systems.")

dt = common_timebase(self.dt, other.dt)

num = polymul(self.num[0][0], other.den[0][0])

den = polymul(self.den[0][0], other.num[0][0])

return TransferFunction(num, den, dt)

# TODO: Division of MIMO transfer function objects is not written yet.

def __rtruediv__(self, other):

"""Right divide two LTI objects."""

if isinstance(other, (int, float, complex, np.number)):

other = _convert_to_transfer_function(

other, inputs=self.ninputs,

outputs=self.ninputs)

else:

other = _convert_to_transfer_function(other)

if (self.ninputs > 1 or self.noutputs > 1 or

other.ninputs > 1 or other.noutputs > 1):

raise NotImplementedError(

"TransferFunction.__rtruediv__ is currently implemented only "

"for SISO systems.")

return other / self

def __pow__(self, other):

if not type(other) == int:

raise ValueError("Exponent must be an integer")

if other == 0:

return TransferFunction([1], [1]) # unity

if other > 0:

return self * (self**(other - 1))

if other < 0:

return (TransferFunction([1], [1]) / self) * (self**(other + 1))

def __getitem__(self, key):

key1, key2 = key

# pre-process

if isinstance(key1, int):

key1 = slice(key1, key1 + 1, 1)

if isinstance(key2, int):

key2 = slice(key2, key2 + 1, 1)

# dim1

start1, stop1, step1 = key1.start, key1.stop, key1.step

if step1 is None:

step1 = 1

if start1 is None:

start1 = 0

if stop1 is None:

stop1 = len(self.num)

# dim1

start2, stop2, step2 = key2.start, key2.stop, key2.step

if step2 is None:

step2 = 1

if start2 is None:

start2 = 0

if stop2 is None:

stop2 = len(self.num[0])

num, den = [], []

for i in range(start1, stop1, step1):

num_i = []

den_i = []

for j in range(start2, stop2, step2):

num_i.append(self.num[i][j])

den_i.append(self.den[i][j])

num.append(num_i)

den.append(den_i)

# Save the label names

outputs = [self.output_labels[i] for i in range(start1, stop1, step1)]

inputs = [self.input_labels[j] for j in range(start2, stop2, step2)]

# Create the system name

sysname = config.defaults['iosys.indexed_system_name_prefix'] + \

self.name + config.defaults['iosys.indexed_system_name_suffix']

return TransferFunction(

num, den, self.dt, inputs=inputs, outputs=outputs, name=sysname)

def freqresp(self, omega):

"""(deprecated) Evaluate transfer function at complex frequencies.

.. deprecated::0.9.0

Method has been given the more pythonic name

:meth:`TransferFunction.frequency_response`. Or use

:func:`freqresp` in the MATLAB compatibility module.

"""

warn("TransferFunction.freqresp(omega) will be removed in a "

"future release of python-control; use "

"sys.frequency_response(omega), or freqresp(sys, omega) in the "

"MATLAB compatibility module instead", DeprecationWarning)

return self.frequency_response(omega)

def poles(self):

"""Compute the poles of a transfer function."""

_, den, denorder = self._common_den(allow_nonproper=True)

rts = []

for d, o in zip(den, denorder):

rts.extend(roots(d[:o + 1]))

return np.array(rts).astype(complex)

def zeros(self):

"""Compute the zeros of a transfer function."""

if self.ninputs > 1 or self.noutputs > 1:

raise NotImplementedError(

"TransferFunction.zeros is currently only implemented "

"for SISO systems.")

else:

# for now, just give zeros of a SISO tf

return roots(self.num[0][0]).astype(complex)

def feedback(self, other=1, sign=-1):

"""Feedback interconnection between two LTI objects."""

other = _convert_to_transfer_function(other)

if (self.ninputs > 1 or self.noutputs > 1 or

other.ninputs > 1 or other.noutputs > 1):

# TODO: MIMO feedback

raise ControlMIMONotImplemented(

"TransferFunction.feedback is currently not implemented for "

"MIMO systems.")

dt = common_timebase(self.dt, other.dt)

num1 = self.num[0][0]

den1 = self.den[0][0]

num2 = other.num[0][0]

den2 = other.den[0][0]

num = polymul(num1, den2)

den = polyadd(polymul(den2, den1), -sign * polymul(num2, num1))

return TransferFunction(num, den, dt)

# For MIMO or SISO systems, the analytic expression is

# self / (1 - sign * other * self)

# But this does not work correctly because the state size will be too

# large.

def minreal(self, tol=None):

"""Remove cancelling pole/zero pairs from a transfer function"""

# based on octave minreal

# default accuracy

from sys import float_info

sqrt_eps = sqrt(float_info.epsilon)

# pre-allocate arrays

num = [[[] for j in range(self.ninputs)] for i in range(self.noutputs)]

den = [[[] for j in range(self.ninputs)] for i in range(self.noutputs)]

for i in range(self.noutputs):

for j in range(self.ninputs):

# split up in zeros, poles and gain

newzeros = []

zeros = roots(self.num[i][j])

poles = roots(self.den[i][j])

gain = self.num[i][j][0] / self.den[i][j][0]

# check all zeros

for z in zeros:

t = tol or \

1000 * max(float_info.epsilon, abs(z) * sqrt_eps)

idx = where(abs(z - poles) < t)[0]

if len(idx):

# cancel this zero against one of the poles

poles = delete(poles, idx[0])

else:

# keep this zero

newzeros.append(z)

# poly([]) returns a scalar, but we always want a 1d array

num[i][j] = np.atleast_1d(gain * real(poly(newzeros)))

den[i][j] = np.atleast_1d(real(poly(poles)))

# end result

return TransferFunction(num, den, self.dt)

def returnScipySignalLTI(self, strict=True):

"""Return a list of a list of :class:`scipy.signal.lti` objects.

For instance,

>>> out = tfobject.returnScipySignalLTI() # doctest: +SKIP

>>> out[3][5] # doctest: +SKIP

is a :class:`scipy.signal.lti` object corresponding to the

transfer function from the 6th input to the 4th output.

Parameters

----------

strict : bool, optional

True (default):

The timebase `tfobject.dt` cannot be None; it must be

continuous (0) or discrete (True or > 0).

False:

if `tfobject.dt` is None, continuous time

:class:`scipy.signal.lti` objects are returned

Returns

-------

out : list of list of :class:`scipy.signal.TransferFunction`

continuous time (inheriting from :class:`scipy.signal.lti`)

or discrete time (inheriting from :class:`scipy.signal.dlti`)

SISO objects

"""

if strict and self.dt is None:

raise ValueError("with strict=True, dt cannot be None")

if self.dt:

kwdt = {'dt': self.dt}

else:

# scipy convention for continuous time lti systems: call without

# dt keyword argument

kwdt = {}

# Preallocate the output.

out = [[[] for j in range(self.ninputs)] for i in range(self.noutputs)]

for i in range(self.noutputs):

for j in range(self.ninputs):

out[i][j] = signalTransferFunction(self.num[i][j],

self.den[i][j],

**kwdt)

return out

def _common_den(self, imag_tol=None, allow_nonproper=False):

"""

Compute MIMO common denominators; return them and adjusted numerators.

This function computes the denominators per input containing all

the poles of sys.den, and reports it as the array den. The

output numerator array num is modified to use the common

denominator for this input/column; the coefficient arrays are also

padded with zeros to be the same size for all num/den.

Parameters

----------

imag_tol: float

Threshold for the imaginary part of a root to use in detecting

complex poles

allow_nonproper : boolean

Do not enforce proper transfer functions

Returns

-------

num: array

n by n by kd where n = max(sys.noutputs,sys.ninputs)

kd = max(denorder)+1

Multi-dimensional array of numerator coefficients. num[i,j]

gives the numerator coefficient array for the ith output and jth

input; padded for use in td04ad ('C' option); matches the

denorder order; highest coefficient starts on the left.

If allow_nonproper=True and the order of a numerator exceeds the

order of the common denominator, num will be returned as None

den: array

sys.ninputs by kd

Multi-dimensional array of coefficients for common denominator

polynomial, one row per input. The array is prepared for use in

slycot td04ad, the first element is the highest-order polynomial

coefficient of s, matching the order in denorder. If denorder <

number of columns in den, the den is padded with zeros.

denorder: array of int, orders of den, one per input

Examples

--------

>>> num, den, denorder = sys._common_den() # doctest: +SKIP

View remainder of file in raw view

FilesExpand file tree

xferfcn.py

Latest commit

History

xferfcn.py

File metadata and controls