Skip to content

[Robust control] Add makeweight function #289

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Closed
wants to merge 2 commits into from
Closed

[Robust control] Add makeweight function #289

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5eb82d1
Fix pep8 warnings
don4get Jan 5, 2019
eb6b0d0
[robust.py] Add method makeweight to compute weighting functions
don4get Apr 13, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
186 changes: 122 additions & 64 deletions control/robust.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,9 +43,13 @@
import numpy as np
from .exception import *
from .statesp import StateSpace
from .xferfcn import tf
from .statefbk import *
import numbers

def h2syn(P,nmeas,ncon):
__all__ = ['h2syn', 'hinfsyn', 'augw', 'mixsyn', 'makeweight']

def h2syn(P, nmeas, ncon):
"""H_2 control synthesis for plant P.

Parameters
Expand All @@ -72,25 +76,25 @@ def h2syn(P,nmeas,ncon):
>>> 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:
# 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)
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]
Expand All @@ -100,7 +104,8 @@ def h2syn(P,nmeas,ncon):

return K

def hinfsyn(P,nmeas,ncon):

def hinfsyn(P, nmeas, ncon):
"""H_{inf} control synthesis for plant P.

Parameters
Expand Down Expand Up @@ -136,26 +141,26 @@ def hinfsyn(P,nmeas,ncon):

"""

#Check for ss system object, need a utility for this?
# 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:
# TODO: Check for continous or discrete, only continuous supported right now
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should raise an exception if a discrete time system is passed as an argument. Use isdtime(P, strict=True) to check.

# 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)
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)
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]
Expand Down Expand Up @@ -202,21 +207,21 @@ def _size_as_needed(w, wname, n):
"""
from . import append, ss
if w is not None:
if not isinstance(w,StateSpace):
if not isinstance(w, StateSpace):
w = ss(w)
if 1==w.inputs and 1==w.outputs:
w = append(*(w,)*n)
if 1 == w.inputs and 1 == w.outputs:
w = append(*(w,) * n)
else:
if w.inputs != n:
msg=("{}: weighting function has {} inputs, expected {}".
format(wname,w.inputs,n))
msg = ("{}: weighting function has {} inputs, expected {}".
format(wname, w.inputs, n))
raise ValueError(msg)
else:
w = ss([],[],[],[])
w = ss([], [], [], [])
return w


def augw(g,w1=None,w2=None,w3=None):
def augw(g, w1=None, w2=None, w3=None):
"""Augment plant for mixed sensitivity problem.

Parameters
Expand All @@ -225,7 +230,6 @@ def augw(g,w1=None,w2=None,w3=None):
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.
Expand Down Expand Up @@ -254,12 +258,12 @@ def augw(g,w1=None,w2=None,w3=None):
ny = g.outputs
nu = g.inputs

w1,w2,w3 = [_size_as_needed(w,wname,n)
for w,wname,n in zip((w1,w2,w3),
('w1','w2','w3'),
(ny,nu,ny))]
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):
if not isinstance(g, StateSpace):
g = ss(g)

# w u
Expand All @@ -270,11 +274,11 @@ def augw(g,w1=None,w2=None,w3=None):
# v [ I | -g ]

# error summer: inputs are -y and r=w
Ie = ss([],[],[],np.eye(ny))
Ie = ss([], [], [], np.eye(ny))
# control: needed to "distribute" control input
Iu = ss([],[],[],np.eye(nu))
Iu = ss([], [], [], np.eye(nu))

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

niw1 = w1.inputs
niw2 = w2.inputs
Expand All @@ -284,43 +288,45 @@ def augw(g,w1=None,w2=None,w3=None):
now2 = w2.outputs
now3 = w3.outputs

q = np.zeros((niw1+niw2+niw3+ny+nu,2))
q[:,0] = np.arange(1,q.shape[0]+1)
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)
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)
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)
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)
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)
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)))
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)
oi = np.arange(1, 1 + now1 + now2 + now3 + ny)

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

return p

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

def mixsyn(g, w1=None, w2=None, w3=None):
"""Mixed-sensitivity H-infinity synthesis.

mixsyn(g,w1,w2,w3) -> k,cl,info
Expand Down Expand Up @@ -356,8 +362,60 @@ def mixsyn(g,w1=None,w2=None,w3=None):
"""
nmeas = g.outputs
ncon = g.inputs
p = augw(g,w1,w2,w3)
p = augw(g, w1, w2, w3)

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


k,cl,gamma,rcond=hinfsyn(p,nmeas,ncon)
info = gamma,rcond
return k,cl,info
def makeweight(dcgain, wc, hfgain, order=1, Ts=None):
""" Compute a weight transfer function, noted W

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need more documentation about the functionality that is implemented. Add a few sentences (and maybe an equation?) describing what the function does.

Also, this function needs to be added to doc/control.rst and doc/matlab.rst.

Parameters
----------
dcgain : double positive scalar
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Incorrect indentation; should line up with Parameters. See numpydocs.

Low frequency gain of 1/W; should be < 1
wc : double strictly positive scalar
Design frequency (where |W| is approximately 1)
hfgain : double positive scalar
High frequency gain of 1/W; should be > 1
order : int strictly positive scalar - optional
Order of the output transfer function
Ts : double strictly positive scalar - optional
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn't the sampling time come from the system (sys.dt)? What happens if they are not compatible?

Sampling time in case W is discrete
W - SISO LTI object
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should be in a section titled Returns. See numpydoc.

"""
if not isinstance(dcgain, numbers.Real):
raise ControlArgument("dcgain should be a real scalar number")
if dcgain < 0:
raise ValueError("dcgain should be a positive real scalar number")

if not isinstance(wc, numbers.Real):
raise ControlArgument("wc should be a real scalar number")
if wc <= 0:
raise ValueError("wc should be a strictly positive real scalar number")

if not isinstance(hfgain, numbers.Real):
raise ControlArgument("hfgain should be a real scalar number")
if hfgain < 0:
raise ValueError("hfgain should be a positive real scalar number")

if not isinstance(order, numbers.Integral):
raise ControlArgument("order should be an integral scalar number")
if order <= 0:
raise ValueError("order should be a strictly positive integral scalar number")

# TODO: Implement discrete behavior
if Ts is not None:
raise ControlNotImplemented("makeweight cannot return discrete weight function")
if Ts is not None and not isinstance(Ts, numbers.Real):
raise ControlArgument("Ts should be a real scalar number")
if Ts is not None and Ts <= 0:
raise ValueError("Ts should be a strictly positive real scalar number")

s = tf([1, 0], [1])

W = (s + wc * (dcgain ** (1. / order))) ** order / (s / (hfgain ** (1. / order)) + wc) ** order

return W
Loading