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 formulationSYS is the discrete-time linear time-invariant model
(A,B,C,D,TS)
TYPE: 2 types,
'reg'
= regulator;'track'
= controller for reference trackingRegulator:
$\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
2x(t+1) = A x(t) + B u(t)
y(t) = C x(t) + D u(t)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
2x(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:
IfCOST.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.
IfCOST.P='lyap'
then P is chosen as the solution of the Lyapunov equationA'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=NcuLIMITS 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:*
IfCOST.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.
IfCOST.P='lyap'
then P is chosen as the solution of the Lyapunov equationA'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 CplexYZEROCON 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-1COST
= cell array of structures of weights:COST{i}.Q
,COST{i}.R
, etc. OnlyCOST{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 constraintsLIMITS{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 | Input: |
refs.y
has as many columns as the number of outputs
1 | Output arguments: |
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 | This syntax is referred to as the 'Single-Output Format'. |
kalman
KALMAN
Design 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 | .xmin, xmax = range of states [ xmin <= x <= xmax ] (linear and 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 | .lpsolver = LP solver (type HELP LPSOL for details) |
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
- Create transfer function model; 2. convert to transfer function model
1 | sys = tf(Numerator,Denominator) |
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 | Numerator = {[1 1] ; 1}; |
Example2
Convert State-Space Model to Transfer Function:
Compute the transfer function of the following state-space model:
1 | sys = ss([-2 -1;1 -2],[1 1;2 -1],[1 0],[0 1]); |
ss
Create state-space model, convert to state-space model
1 | sys = ss(A,B,C,D) |
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 anNy
-by-Nx
real- or complex-valued matrix.D
is anNy
-by-Nu
real- or complex-valued matrix.To set
D = 0
, setD
to the scalar0
(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 | sysd = c2d(sys,Ts) |
step
Step response plot of dynamic system; step response data
1 | step(sys) |
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 function1
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]
minimizes1
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 | K = 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 A – BK 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 A – BK 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 | MPCOBJ = mpc(MODELS,TS,P,M) where MODELS is a structure containing |
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 | f = msgbox(message) |
warndlg
Create warning dialog box
1 | f = warndlg(msg) |
waitfor
waitfor
: Block execution and wait for condition [after operation, will terminate the message box]
1 | waitfor(mo) |
Examples:
1 | mydlg = warndlg('This is a warning.', 'A Warning Dialog'); |
1 | gcf:get current figure;%当前figure |