# robust.py - tools for robust control

#

# Author: Steve Brunton, Kevin Chen, Lauren Padilla

# Date: 24 Dec 2010

#

# This file contains routines for obtaining reduced order models

#

# Copyright (c) 2010 by California Institute of Technology

# All rights reserved.

#

# 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.

#

# $Id$

# External packages and modules

import numpy as np

from .exception import *

from .statesp import StateSpace

from .statefbk import *

def h2syn(P, nmeas, ncon):

"""H_2 control synthesis for plant P.

Parameters

----------

P: partitioned lti plant (State-space sys)

nmeas: number of measurements (input to controller)

ncon: number of control inputs (output from controller)

Returns

-------

K: controller to stabilize P (State-space sys)

Raises

------

ImportError

if slycot routine sb10hd is not loaded

See Also

--------

StateSpace

Examples

--------

>>> K = h2syn(P,nmeas,ncon)

"""

# Check for ss system object, need a utility for this?

# TODO: Check for continous or discrete, only continuous supported right now

# if isCont():

# dico = 'C'

# elif isDisc():

# dico = 'D'

# else:

dico = 'C'

try:

from slycot import sb10hd

except ImportError:

raise ControlSlycot("can't find slycot subroutine sb10hd")

n = np.size(P.A, 0)

m = np.size(P.B, 1)

np_ = np.size(P.C, 0)

out = sb10hd(n, m, np_, ncon, nmeas, P.A, P.B, P.C, P.D)

Ak = out[0]

Bk = out[1]

Ck = out[2]

Dk = out[3]

K = StateSpace(Ak, Bk, Ck, Dk)

return K

def hinfsyn(P, nmeas, ncon):

"""H_{inf} control synthesis for plant P.

Parameters

----------

P: partitioned lti plant

nmeas: number of measurements (input to controller)

ncon: number of control inputs (output from controller)

Returns

-------

K: controller to stabilize P (State-space sys)

CL: closed loop system (State-space sys)

gam: infinity norm of closed loop system

rcond: 4-vector, reciprocal condition estimates of:

1: control transformation matrix

2: measurement transformation matrix

3: X-Riccati equation

4: Y-Riccati equation

TODO: document significance of rcond

Raises

------

ImportError

if slycot routine sb10ad is not loaded

See Also

--------

StateSpace

Examples

--------

>>> K, CL, gam, rcond = hinfsyn(P,nmeas,ncon)

"""

# Check for ss system object, need a utility for this?

# TODO: Check for continous or discrete, only continuous supported right now

# if isCont():

# dico = 'C'

# elif isDisc():

# dico = 'D'

# else:

dico = 'C'

try:

from slycot import sb10ad

except ImportError:

raise ControlSlycot("can't find slycot subroutine sb10ad")

n = np.size(P.A, 0)

m = np.size(P.B, 1)

np_ = np.size(P.C, 0)

gamma = 1.e100

out = sb10ad(n, m, np_, ncon, nmeas, gamma, P.A, P.B, P.C, P.D)

gam = out[0]

Ak = out[1]

Bk = out[2]

Ck = out[3]

Dk = out[4]

Ac = out[5]

Bc = out[6]

Cc = out[7]

Dc = out[8]

rcond = out[9]

K = StateSpace(Ak, Bk, Ck, Dk)

CL = StateSpace(Ac, Bc, Cc, Dc)

return K, CL, gam, rcond

def _size_as_needed(w, wname, n):

"""Convert LTI object to appropriately sized StateSpace object.

Intended for use in .robust only

Parameters

----------

w: None, 1x1 LTI object, or mxn LTI object

wname: name of w, for error message

n: number of inputs to w

Returns

-------

w_: processed weighting function, a StateSpace object:

- if w is None, empty StateSpace object

- if w is scalar, w_ will be w * eye(n)

- otherwise, w as StateSpace object

Raises

------

ValueError

- if w is not None or scalar, and doesn't have n inputs

See Also

--------

augw

"""

from . import append, ss

if w is not None:

if not isinstance(w, StateSpace):

w = ss(w)

if 1 == w.ninputs and 1 == w.noutputs:

w = append(*(w,) * n)

else:

if w.ninputs != n:

msg = ("{}: weighting function has {} inputs, expected {}".

format(wname, w.ninputs, n))

raise ValueError(msg)

else:

w = ss([], [], [], [])

return w

def augw(g, w1=None, w2=None, w3=None):

"""Augment plant for mixed sensitivity problem.

Parameters

----------

g: LTI object, ny-by-nu

w1: weighting on S; None, scalar, or k1-by-ny LTI object

w2: weighting on KS; None, scalar, or k2-by-nu LTI object

w3: weighting on T; None, scalar, or k3-by-ny LTI object

p: augmented plant; StateSpace object

If a weighting is None, no augmentation is done for it. At least

one weighting must not be None.

If a weighting w is scalar, it will be replaced by I*w, where I is

ny-by-ny for w1 and w3, and nu-by-nu for w2.

Returns

-------

p: plant augmented with weightings, suitable for submission to hinfsyn or h2syn.

Raises

------

ValueError

- if all weightings are None

See Also

--------

h2syn, hinfsyn, mixsyn

"""

from . import append, ss, connect

if w1 is None and w2 is None and w3 is None:

raise ValueError("At least one weighting must not be None")

ny = g.noutputs

nu = g.ninputs

w1, w2, w3 = [_size_as_needed(w, wname, n)

for w, wname, n in zip((w1, w2, w3),

('w1', 'w2', 'w3'),

(ny, nu, ny))]

if not isinstance(g, StateSpace):

g = ss(g)

# w u

# z1 [ w1 | -w1*g ]

# z2 [ 0 | w2 ]

# z3 [ 0 | w3*g ]

# [------+--------- ]

# v [ I | -g ]

# error summer: inputs are -y and r=w

Ie = ss([], [], [], np.eye(ny))

# control: needed to "distribute" control input

Iu = ss([], [], [], np.eye(nu))

sysall = append(w1, w2, w3, Ie, g, Iu)

niw1 = w1.ninputs

niw2 = w2.ninputs

niw3 = w3.ninputs

now1 = w1.noutputs

now2 = w2.noutputs

now3 = w3.noutputs

q = np.zeros((niw1 + niw2 + niw3 + ny + nu, 2))

q[:, 0] = np.arange(1, q.shape[0] + 1)

# Ie -> w1

q[:niw1, 1] = np.arange(1 + now1 + now2 + now3,

1 + now1 + now2 + now3 + niw1)

# Iu -> w2

q[niw1:niw1 + niw2, 1] = np.arange(1 + now1 + now2 + now3 + 2 * ny,

1 + now1 + now2 + now3 + 2 * ny + niw2)

# y -> w3

q[niw1 + niw2:niw1 + niw2 + niw3, 1] = np.arange(1 + now1 + now2 + now3 + ny,

1 + now1 + now2 + now3 + ny + niw3)

# -y -> Iy; note the leading -

q[niw1 + niw2 + niw3:niw1 + niw2 + niw3 + ny, 1] = -np.arange(1 + now1 + now2 + now3 + ny,

1 + now1 + now2 + now3 + 2 * ny)

# Iu -> G

q[niw1 + niw2 + niw3 + ny:niw1 + niw2 + niw3 + ny + nu, 1] = np.arange(

1 + now1 + now2 + now3 + 2 * ny,

1 + now1 + now2 + now3 + 2 * ny + nu)

# input indices: to Ie and Iu

ii = np.hstack((np.arange(1 + now1 + now2 + now3,

1 + now1 + now2 + now3 + ny),

np.arange(1 + now1 + now2 + now3 + ny + nu,

1 + now1 + now2 + now3 + ny + nu + nu)))

# output indices

oi = np.arange(1, 1 + now1 + now2 + now3 + ny)

p = connect(sysall, q, ii, oi)

return p

def mixsyn(g, w1=None, w2=None, w3=None):

"""Mixed-sensitivity H-infinity synthesis.

mixsyn(g,w1,w2,w3) -> k,cl,info

Parameters

----------

g: LTI; the plant for which controller must be synthesized

w1: weighting on s = (1+g*k)**-1; None, or scalar or k1-by-ny LTI

w2: weighting on k*s; None, or scalar or k2-by-nu LTI

w3: weighting on t = g*k*(1+g*k)**-1; None, or scalar or k3-by-ny LTI

At least one of w1, w2, and w3 must not be None.

Returns

-------

k: synthesized controller; StateSpace object

cl: closed system mapping evaluation inputs to evaluation outputs; if

p is the augmented plant, with

[z] = [p11 p12] [w],

[y] [p21 g] [u]

then cl is the system from w->z with u=-k*y. StateSpace object.

info: tuple with entries, in order,

- gamma: scalar; H-infinity norm of cl

- rcond: array; estimates of reciprocal condition numbers

computed during synthesis. See hinfsyn for details

If a weighting w is scalar, it will be replaced by I*w, where I is

ny-by-ny for w1 and w3, and nu-by-nu for w2.

See Also

--------

hinfsyn, augw

"""

nmeas = g.noutputs

ncon = g.ninputs

p = augw(g, w1, w2, w3)

k, cl, gamma, rcond = hinfsyn(p, nmeas, ncon)

info = gamma, rcond

return k, cl, info