Hybrid Toolbox

linear system

###lincon

lincon Constrained linear quadratic controller for linear systems

  • L=lincon(SYS,TYPE,COST,INTERVAL,LIMITS,QPSOLVER) build a regulator (or a controller for reference tracking) based on a finite time linear quadratic constrained optimal control formulation

    • SYS is the discrete-time linear time-invariant model (A,B,C,D,TS)

    • TYPE: 2 types, 'reg' = regulator; 'track' = controller for reference tracking

      1. Regulator:

        $\min x’(N)Px(N) + \sum\limits_{k=0}^{N-1} x’(k)Qx(k) + u’(k)Ru(k) + \rho*\epsilon ^2$

        s.t. $ y_{min} - \epsilon <= y(k) <=y_{max} + \epsilon, k=1,…,N_{cy}$
        ​ $u_{min} <= u(k) <=u_{max}, k=0,…,N_{cu}$
        ​ $u(k) = Kx(k), k>=N_u$

      1
      2
      x(t+1) = A x(t) + B u(t)
      y(t) = C x(t) + D u(t)
      1. Reference Tracking :

        $\min \sum\limits_{k=0}^{N-1} [y(k)-r]’S[y(k)-r] + \Delta u’(k) T \Delta u(k) + \rho*\epsilon^2$

        s.t. $y_{min} - \epsilon <= y(k) <=y_{max} + \epsilon,\ k=1,…,N_{cy}$
        ​ $u_{min} <= u(k) <=u_{max}, \ k=0,…,N_{cu}$
        ​ $\Delta u_{min} <= du(k) <=\Delta u_{max}, \ k=0,…,N_{cu}$
        ​ $\Delta u(k) = 0, \ k>=N_u$

      1
      2
      x(t+1) = Ax(t) + B [u(t-1)+du(t)]   %du(t) = u(t)-u(t-1) = input increment
      y(t) = Cx(t) + D [u(t-1)+du(t)]
    • COST is a structure of weights with fields:

      1
      2
      3
      4
      5
      6
      7
      .Q = weight Q on states             [ x'(k)Qx(k) ]           (regulator only)
      .R = weight R on inputs [ u'(k)Ru(k) ] (regulator only)
      .P = weight P on final state [ x'(N)Px(N) ] (regulator only)
      .S = weight S on outputs [ (y(k)-r)'S(y(k)-r) ] (tracking only)
      .T = weight T on input increments [ du'(k) T du(k) ] (tracking only)
      .rho = weight on slack var [ rho*eps^2 ] (regulator and tracking) %rho=Inf means hard output constraints
      .K = feedback gain (default: K=0) [ u(k)=K*x(k) for k>=Nu] (only regulator)

      Only for regulators:
      If COST.P='lqr' then P is chosen as the solution of the LQR problem with weights Q,R, and K is chosen as the corresponding LQR gain.
      If COST.P='lyap' then P is chosen as the solution of the Lyapunov equation A'PA - P + Q = 0, and K=0 (default).

    • INTERVAL is a structure of number of input and output optimal control steps

      1
      2
      3
      4
      .N   = optimal control interval over which the cost function is summed
      .Nu = number of free optimal control moves u(0),...,u(Nu-1)
      .Ncy = output constraints are checked up to time k=Ncy
      .Ncu = input constraints are checked up to time k=Ncu
    • LIMITS is a structure of constraints with fields:

      1
      2
      3
      4
      5
      6
      .umin = lower bounds on inputs             [ u(k)>=umin ]
      .umax = upper bounds on inputs [ u(k)<=umax ]
      .ymin = lower bounds on outputs [ y(k)>=ymin ]
      .ymax = upper bounds on outputs [ y(k)<=ymax ]
      .dumin = lower bounds on input increments [ du(k)>=dumin ] (only tracking)
      .dumax = upper bounds on input increments [ du(k)<=dumax ] (only tracking)

      Only for regulators:*
      If COST.P='lqr' then P is chosen as the solution of the LQR problem with weights Q,R, and K is chosen as the corresponding LQR gain.
      If COST.P='lyap' then P is chosen as the solution of the Lyapunov equation A'PA - P + Q = 0, and K=0 (default).

    • QPSOLVER specifies the QP used for evaluating the control law. Valid options are:

      1
      2
      3
      4
      5
      'qpact'      active set method
      'quadprog' QUADPROG from Optimization Toolbox
      'qp' QP from Optimization Toolbox
      'nag' QP from NAG Foundation Toolbox
      'cplex' QP from Ilog Cplex
    • YZEROCON enforce output constraints also at prediction time k=0 if YZEROCON=1 (default: YZEROCON=0).

      If YEZEROCON is a 0/1 vector of the same dimension as the output vector, then only those outputs where YZEROCON(i)=1 are constrained at time k=0.

    Time-varying formulation: Use the same syntax as above, with input arguments as follow:

  • SYS = cell array of linear models, SYS{i}=model used to predict at step t+i-1

  • COST = cell array of structures of weights: COST{i}.Q, COST{i}.R, etc. Only COST{1}.rho is used for penalty soft constraint penalty, COST{i}.rho is ignored for i>1.

  • LIMITS = cell arrays of structures of upper and lower bounds: LIMITS{i}.umin, etc. The constraints LIMITS{i} refer to time index t+i, for both inputs and outputs.linsfun

sim for lincon

sim : Closed-loop simulation of constrained optimal controllers for linear systems.

[X,U,T,Y]=sim(LINCON,MODEL,refs,x0,Tstop,u0,qpsolver,verbose) simulates the closed-loop of the linear system MODEL with the controller LINCON based on quadratic programming.

Usually, LINCON is based on model MODEL (nominal closed-loop).

1
2
3
4
5
6
7
8
9
Input: 

LINCON = optimal controller
refs = structure or references with fields; ref.y output references (only for tracking)
x0 = initial state
Tstop = total simulation time (total number of steps = Tstop/LINCON.Ts)
u0 = input at time t=-1 (only for tracking)
qpsolver = QP solver (type HELP QPSOL for details)
verbose = verbosity level (0,1)

refs.y has as many columns as the number of outputs

1
2
3
4
5
6
Output arguments:
X is the sequence of states, with as many columns as the number of states.
U is the sequence of inputs, with as many columns as the number of inputs.
T is the vector of time steps.
I is the sequence of region numbers
Y is the sequence of outputs, with as many columns as the number of outputs.

Without output arguments, sim produces a plot of the trajectories

LINCON objects based on time-varying prediction models/weights/limits can be simulated in closed-loop with LTI models. Time-varying simulation models are not supported, use custom for-loops instead.

sim (general)

sim : Simulate a Simulink model

SimOut = sim('MODEL', PARAMETERS) simulates your Simulink model,

'PARAMETERS' represents a list of parameter name-value pairs, a structure containing parameter settings, or a configuration set.

SimOut returned by the sim command is an object that contains all of the logged simulation results.

PARAMETERS can be used to override existing block diagram configuration parameters for the duration of the simulation.

1
2
3
4
5
6
7
8
This syntax is referred to as the 'Single-Output Format'.
paramNameValStruct.SimulationMode = 'rapid';
paramNameValStruct.AbsTol = '1e-5';
paramNameValStruct.SaveState = 'on';
paramNameValStruct.StateSaveName = 'xoutNew';
paramNameValStruct.SaveOutput = 'on';
paramNameValStruct.OutputSaveName = 'youtNew';
simOut = sim('vdp',paramNameValStruct)

kalman

KALMANDesign Kalman filter for constrained optimal controllers

KALMAN(CON,Q,R) design a state observer for the linear model (A,B,C,D)=CON.model which controller CON is based on, using Kalman filtering techniques. The resulting observer is stored in the 'Observer' field of the controller object. The controller object can be either of class LINCON or EXPCON.

Q is the covariance matrix of state noise, R is the covariance matrix of outputnoise.

Kest=KALMAN(CON,Q,R) also return the state observer as an LTI object.

Kest receives u[n] and y[n] as inputs (in this order), and provides the best estimated state x[n|n-1] of x[n] given the past measurements y[n-1], y[n-2],...

1
x[n+1|n] = Ax[n|n-1] + Bu[n] + L(y[n] - Cx[n|n-1] - Du[n])

Kest=KALMAN(CON,Q,R,ymeasured) assumes that only the outputs specified in vector ymeasured are measurable outputs. In this case, R is a nym-by-nym matrix, where nym is the number of measured outputs.

[Kest,M]=KALMAN(CON,Q,R,ymeasured) also returns the gain M which allows computing the measurement update

1
x[n|n]  = x[n|n-1] + M(y[n] - Cx[n|n-1] - Du[n])

For time-varying prediction models, the first model CON.model{1} is used to build the Kalman filter, as it is the model of the process for the current time step.

expcon

EXPCON Explicit controllers (for linear/hybrid systems)

**E=EXPCON(C,RANGE,OPTIONS) **convert controller C to piecewise affine explicit form.

C is a constrained optimal controller for linear systems or hybrid systems (an object of class @LINCON or @HYBCON)

RANGE defines the range of the parameters the explicit solution is found with respect to (it may be different from the limits over variables imposed in the control law). RANGE is a structure with fields:

1
2
3
4
5
.xmin, xmax = range of states            [ xmin <= x <= xmax ]       (linear and hybrid)
.umin, umax = range of inputs [ umin <= u <= umax ] (linear w/ tracking)
.refymin, refymax = range of output refs [ refymin <= r.y <= refymax] (linear w/ tracking and hybrid)
.refxmin, refxmax = range of state refs [ refxmin <= r.x <= refxmax ] (hybrid)
.refumin, refumax = range of input refs [ refumin <= r.u <= refumax ] (hybrid)

For hybrid controllers, refxmin,refxmax,refumin,refumax,refymin,refymax must have dimensions consistent with the number of indices specified in C.ref (type “help hybcon” for details)

OPTIONS defines various options for computing the explicit control law:

1
2
3
4
5
6
7
8
9
10
11
12
.lpsolver    = LP solver (type HELP LPSOL for details)
.qpsolver = QP solver (type HELP QPSOL for details)
.fixref = structure with fields 'y','x','u' defining references that are fixed at given values, therefore simplifying the parametric solution [hybrid only]
.valueref = structure with fields 'y','x','u' of values at which references are fixed [hybrid only]
.flattol Tolerance for a polyhedral set to be considered flat
.waitbar display waitbar (only for hybrid)
.verbose = level of verbosity {0,1,2} of mp-PWA solver
.mplpverbose = level of verbosity {0,1,2} of mp-LP solver
.uniteeps = tolerance for judging convexity of the union of polyhedra
.join = flag for reducing the complexity by joining regions whose union is a convex set.
.reltol = tolerance used for several polyhedral operations
.sequence = flag for keeping the entire optimal control sequence as an explicit function of the parameters also in the linear case. Flag ignored in the hybrid case: sequence always stored with quadratic costs, never stored with inf-norms.

Example: to fix the reference signal for x(1),x(2) at the values rx(1)=0.6, rx(2)=-1.4 and mantain the reference rx(3) for x(3) as a free parameter, specify options.fixref.x=[1 2], options.valueref.x=[0.6 -1.4].

Note: @MPC objects from MPC Toolbox are no longer supported, as they are handled directly in the MPC Toolbox since R2014b.

See also EXPCON/EVAL for the order states and references are stored inside the parameter vector.

##Modelling

tf

  1. Create transfer function model; 2. convert to transfer function model
1
2
3
4
5
6
7
8
sys = tf(Numerator,Denominator) 
sys = tf(Numerator,Denominator,Ts)
sys = tf(M)
sys = tf(Numerator,Denominator,ltisys)
tfsys = tf(sys)
tfsys = tf(sys, 'measured')
tfsys = tf(sys, 'noise')
tfsys = tf(sys, 'augmented')

Example1:

Create Transfer Function with One Input and Two Outputs:

​ Create the following transfer function model:

​ $$H(p)=\begin{bmatrix} \frac{p+1}{p^2+2*p+2} \ \frac{1}{p}\end{bmatrix}$$

1
2
3
4
5
Numerator = {[1 1] ; 1};
Denominator = {[1 2 2] ; [1 0]};
H = tf(Numerator,Denominator,'InputName','current',...
'OutputName',{'torque' 'ang. velocity'},...
'Variable','p')

Example2

Convert State-Space Model to Transfer Function:

​ Compute the transfer function of the following state-space model:

1
2
sys = ss([-2 -1;1 -2],[1 1;2 -1],[1 0],[0 1]);
tf(sys)

ss

Create state-space model, convert to state-space model

1
2
3
4
5
6
7
8
9
10
sys = ss(A,B,C,D) 
sys = ss(A,B,C,D,Ts)
sys = ss(D)
sys = ss(A,B,C,D,ltisys)
sys_ss = ss(sys)
sys_ss = ss(sys,'minimal')
sys_ss = ss(sys,'explicit')
sys_ss = ss(sys, 'measured')
sys_ss = ss(sys, 'noise')
sys_ss = ss(sys, 'augmented')

sys = ss(A,B,C,D) creates a state-space model object representing the continuous-time state-space model

​ $x=Ax+Bu$

​ $y=Cx+Du$

​ For a model with Nx states, Ny outputs, and Nu inputs:

​ - A is an Nx-by-Nx real- or complex-valued matrix.

​ - B is an Nx-by-Nu real- or complex-valued matrix.

  • C is an Ny-by-Nx real- or complex-valued matrix.

  • D is an Ny-by-Nu real- or complex-valued matrix.

    To set D = 0 , set D to the scalar 0 (zero), regardless of the dimension.

sys = ss(A,B,C,D,Ts) creates the discrete-time model

​ $x[n+1]=Ax[n]+Bu[n]$

​ $y[n]=Cx[n]+Du[n]$

with sample time Ts (in seconds). Set Ts = -1 or Ts = [] to leave the sample time unspecified.

ssdata

ssdata : Quick access to state-space data.

[A,B,C,D] = ssdata(SYS) returns the A,B,C,D matrices of the state-space model SYS. If SYS is not a state-space model, it is first converted to the state-space representation. If SYS is in descriptor form (nonempty
E matrix), an equivalent explicit form is first derived. If SYS has internal delays, A,B,C,D are obtained by first setting all internal delays to zero (delay-free dynamics).

[A,B,C,D,TS] = ssdata(SYS) also returns the sample time TS. Other properties of SYS can be accessed using struct-like dot syntax (for example, SYS.StateName).

Support for arrays of state-space models:

[A,B,C,D,TS] = ssdata(SYS,J1,...,JN) extracts the data for the (J1,…,JN) entry in the model array SYS.

If all models in SYS have the same order, ssdata returns multi-dimensional arrays A,B,C,D where A(:,:,k), B(:,:,k), C(:,:,k), D(:,:,k) are the state-space matrices of the k-th model SYS(:,:,k).

If the models in SYS have variable order, use the syntax [A,B,C,D] = ssdata(SYS,'cell') to extract the state-space matrices of each model as separate cells in the cell arrays A,B,C,D.

c2d

Convert model from continuous to discrete time

1
2
3
4
5
sysd = c2d(sys,Ts)
sysd = c2d(sys,Ts,method)
sysd = c2d(sys,Ts,opts)
[sysd,G] = c2d(sys,Ts,method)
[sysd,G] = c2d(sys,Ts,opts)

step

Step response plot of dynamic system; step response data

1
2
3
4
5
6
7
8
9
10
11
12
step(sys) 
step(sys,Tfinal)
step(sys,t)
step(sys1,sys2,...,sysN)
step(sys1,sys2,...,sysN,Tfinal)
step(sys1,sys2,...,sysN,t)
y = step(sys,t)
[y,t] = step(sys)
[y,t] = step(sys,Tfinal)
[y,t,x] = step(sys)
[y,t,x,ysd] = step(sys)
[y,...] = step(sys,...,options)

solution

dare

dare : Solve discrete-time algebraic Riccati equations.

[X,L,G] = dare(A,B,Q,R,S,E) computes the unique stabilizing solution X of the discrete-time algebraic Riccati equation $E’XE = A’XA - (A’XB + S)(B’XB + R)^{-1} (A’XB + S)’ + Q$ or, equivalently (if R is nonsingular) $E’XE = F’XF - F’XB(B’XB + R)^{-1} B’XF + Q - SR^{-1} S’$ with $F:=A-BR^{-1}S’$.

When omitted, $R$, $S$ and $E$ are set to the default values $R=I$ , $S=0$ , and $E=I$. Beside the solution X, dare also returns the gain matrix $G = (B’XB + R)^{-1}(B’XA + S’)$, and the vector $L$ of closed-loop eigenvalues (i.e., EIG(A-B*G,E)).

[X,L,G,REPORT] = dare(...) returns a diagnosis REPORT with value:

  • -1 if the Symplectic matrix has eigenvalues on the unit circle
  • -2 if there is no finite stabilizing solution X
  • Frobenius norm of relative residual if X exists and is finite. This syntax does not issue any error message when X fails to exist.

[X1,X2,D,L] = dare(A,B,Q,...,'factor') returns two matrices X1, X2 and a diagonal scaling matrix D such that $X = D(X_2/X_1)D$. The vector $L$ contains the closed-loop eigenvalues. All outputs are empty when the associated Symplectic matrix has eigenvalues on the unit circle

dlqr

dlqr : Linear-quadratic regulator design for discrete-time systems.

[K,S,E] = dlqr(A,B,Q,R,N) calculates the optimal gain matrix $K$ such that the state-feedback law u[n] = -Kx[n] minimizes the cost function

  J = Sum {x'Qx + u'Ru + 2*x'Nu}

subject to the state dynamics x[n+1] = Ax[n] + Bu[n].

The matrix N is set to zero when omitted. Also returned are the Riccati equation solution S and the closed-loop eigenvalues E :

1
A'SA - S - (A'SB+N) (R+B'SB)^(-1) (B'SA+N') + Q = 0,   E = EIG(A-B*K)

lqr

lqr : Linear-quadratic regulator design for state space systems.

[K,S,E] = lqr(SYS,Q,R,N) calculates the optimal gain matrix K such that:

  • For a continuous-time state-space model SYS, the state-feedback law u = -Kx minimizes the cost function

    1
    J = Integral {x'Qx + u'Ru + 2*x'Nu} dt

    subject to the system dynamics dx/dt = Ax + Bu

  • For a discrete-time state-space model SYS, u[n] = -Kx[n] minimizes

    1
    J = Sum {x'Qx + u'Ru + 2*x'Nu}

    subject to x[n+1] = Ax[n] + Bu[n].

The matrix N is set to zero when omitted. Also returned are the solution S of the associated algebraic Riccati equation and the closed-loop eigenvalues E = EIG(A-B*K) .

[K,S,E] = lqr(A,B,Q,R,N) is an equivalent syntax for continuous-time models with dynamics dx/dt = Ax + Bu

operspec

operspec : Create operating point specifications for Simulink model.

OP = operspec('sys') returns an operating point specification object OP, for a Simulink model 'sys'. Edit the default operating point specifications directly or use get and set. Use the operating point specifications object with the function FINDOP to find operating points based on the specifications. Use these operating points with the function LINEARIZE to create linearized models.

OP = operspec('sys', N) returns an N-dimensional column vector of the default operating point specification objects for a Simulink model 'sys'.

OP = operspec('sys', SIZE) returns an array of the default operating point specification objects. The size of OP is specified in the SIZE vector. For example, operspec('vdp', [2,3]) returns a 2-by-3 matrix of default operating point specifications for the Simulink model 'vdp'.

findop

findop : Find operating points from specifications or simulation

place

place Pole placement design

1
2
K = place(A,B,p)
[K,prec,message] = place(A,B,p)

Given the single- or multi-input system $x=Ax+Bu$, and a vector p of desired self-conjugate closed-loop pole locations, place computes a gain matrix K such that the state feedback u=-Kx places the closed-loop poles at the locations p. In other words, the eigenvalues of ABK match the entries of p (up to the ordering).

K = place(A,B,p) places the desired closed-loop poles p by computing a state-feedback gain matrix K. All the inputs of the plant are assumed to be control inputs. The length of p must match the row size of A. place works for multi-input systems and is based on the algorithm from [1]. This algorithm uses the extra degrees of freedom to find a solution that minimizes the sensitivity of the closed-loop poles to perturbations in A or B.

[K,prec,message] = place(A,B,p) returns prec, an estimate of how closely the eigenvalues of ABK match the specified locations p (prec measures the number of accurate decimal digits in the actual closed-loop poles). If some nonzero closed-loop pole is more than 10% off from the desired location, message contains a warning message.

mpc

mpc Create an mpc controller.

MPCOBJ = mpc(PLANT) creates the mpc object based on the linear plant model PLANT, which must be a SS, TF, ZPK or a linear model from System Identification Toolbox such as IDSS, IDTF, IDPROC, IDPOLY and IDGREY.

MPCOBJ = mpc(PLANT,TS) specifies a sampling time for the mpc controller. A continuous-time PLANT model is discretized with sampling time TS. A discrete-time PLANT model is resampled to sampling time TS.

MPCOBJ = mpc(PLANT,TS,P) specifies the prediction horizon P. P is a positive, finite integer.

MPCOBJ = mpc(PLANT,TS,P,M) specifies the control horizon M. M is either an integer (<= P) or a vector of blocked moves such that sum(M)<= P.

MPCOBJ = mpc(PLANT,TS,P,M,Weights,MV,OV,DV) also allows you to specify Weights, constraints on MV, MV rates and plant outputs, scale factors, as well as other properties with structures Weights, MV, OV and DV in that order. Those properties can also be specified using either

1
SET(MPCOBJ,Property1,Value1,Property2,Value2,...)

or

1
MPCOBJ.Property = Value

Disturbance models and operating conditions can be specified using the syntax:

1
2
3
4
5
6
7
MPCOBJ = mpc(MODELS,TS,P,M) where MODELS is a structure containing 
MODELS.Plant = plant model (LTI or linear model from System Identification Toolbox).
.Disturbance = model describing the input disturbances. This model is assumed to be driven by unit variance white noise.
.Noise = model describing the plant output measurement noise. This model is assumed to be driven by unit variance white noise.
.Nominal = structure specifying the plant nominal conditions at which the plant was linearized.

MPCOBJ = mpc creates an empty mpc object

initial

initial Initial condition response of state-space models.

initial(SYS,X0) plots the undriven response of the state-space model SYS (created with SS) with initial condition X0 on the states. This response is characterized by the equations
.
Continuous time: $ x = A x ,\ y = C x ,\ x(0) = x_0 $

Discrete time: $x(k+1) = Ax(k), \ y(k) = Cx(k), \ x(0) = x_0 $

The time range and number of points are chosen automatically.

initial(SYS,X0,T) uses the time vector T for simulation (expressed in the time units of SYS). For discrete-time models, T should be of the form Ti:Ts:Tf where Ts is the sample time. For continuous-time models, T should be of the form Ti:dt:Tf where dt is the sampling period for the discrete approximation of SYS.

initial(SYS,X0,T) uses the time vector T for simulation (expressed in the time units of SYS). For discrete-time models, T should be of the form Ti:Ts:Tf where Ts is the sample time. For continuous-time models, T should be of the form Ti:dt:Tf where dt is the sampling period for the discrete approximation of SYS.

run sim

lsim

lsim Simulate time response of dynamic systems to arbitrary inputs.

lsim(SYS,U,T) plots the time response of the dynamic system SYS to the input signal described by U and T. The time vector T is expressed in the time units of SYS and consists of regularly spaced time samples. The matrix U has as many columns as inputs in SYS and its i-th row specifies the input value at time T(i). For example, t = 0:0.01:5; u = sin(t); lsim(sys,u,t) simulates the response of a single-input model SYS to the input u(t)=sin(t) during 5 time units.

For discrete-time models, U should be sampled at the same rate as SYS (T is then redundant and can be omitted or set to the empty matrix). For continuous-time models, choose the sampling period T(2)-T(1) small enough to accurately describe the input U. lsim issues a warning when U is undersampled and hidden oscillations may occur.

msgbox

msgbox : Create message dialog box

1
2
3
4
5
f = msgbox(message)
f = msgbox(message,title)
f = msgbox(message,title,icon)
f = msgbox(message,title,'custom',icondata,iconcmap)
f = msgbox(___,createmode)

warndlg

Create warning dialog box

1
2
3
4
f = warndlg(msg)
f = warndlg(msg,title)
f = warndlg(msg,title,opts)
f = warndlg

waitfor

waitfor : Block execution and wait for condition [after operation, will terminate the message box]

1
2
3
waitfor(mo)
waitfor(mo,propname)
waitfor(mo,propname,propvalue)

Examples:

1
2
3
mydlg = warndlg('This is a warning.', 'A Warning Dialog');
waitfor(mydlg);
disp('This prints after you close the warning dialog.');
1
2
3
gcf:get current figure%当前figure
gca:get current axes;%当前坐标轴
gcbo:get callback object %调用callback的对象句柄