Models of physical processes Lecture notes 5DD20 January, 2006
P.P.J. van den Bosch EH4.33
[email protected]
Contents Processes and Models
1. Modelling 1.1 Models 1.1.1 Differential equations 1.1.2 State space equations 1.2 Deriving a model 1.2.1 SI units 1.2.2 Analogies among different physical fields 1.2.3 Deriving a model in 4 steps 1.2.4 Sign rule 1.3 Models in different physical fields 1.3.1 Electrical processes 1.3.2 Mechanical processes 1.3.3 Hydraulic processes 1.3.4 Thermal processes 1.3.5 Mass flow, compartment models (pharmacy kinetics) 1.3.6 Chemical reactions 1.4 Linearisation of non-linear state equations
2. Transfer functions 2.1 Laplace Transformation 2.2 Relations between differential equation, state space model and transfer function 2.2.1 Derivation of a state space model from a differential equation 2.2.2 Derivation of the transfer function from a state space model 2.2.3 Derivation of the differential equation from a transfer function 2.3 Frequency response
3. Regulation 3.1 Static feedback loop 3.2 Dynamic feedback loop 3.3 Steady state error
4. Exercises 1
Processes and models These condensed lecture notes deal with models of physical, for example electrical, mechanical, hydraulic, thermal, pharmakinetic and chemical processes. It turns out that these different types of processes cab ne represented by the same kind of dynamical (mathematical) model can be derived. By using differential equations and state space models the same mathematical formulations are used. Deriving these mathematical models from a physical process is called modelling, as illustrated in figure 1. For analysing the dynamic behaviour of these mathematical models and, consequently, of the underlying physical processes, the transfer function is introduced. It turns out that the characteristics of these transfer functions can be interpreted directly to predict its dynamic behaviour, for example its step response.
u
Process
Model
y(t)
y
t
Modelling
Analysis
Figure 1. Modelling and analysis
2
1. Modelling 1.1 Models Models are used to understand and predict the behaviour of real-life processes. As soon as we can describe these real-life processes sufficiently accurate with the aid of a mathematical model, analysing and understanding these processes becomes increasingly easier. When the model is even sufficiently accurate, it is possible to make predictions on future behaviour of the process (e.g. weather forecasting). As such, models can be considered as our knowledge of real-life processes. One of the most difficult steps is to determine the system boundary in the real-life process with inside the components and functions to be modelled and outside other parts that have to be neglected. If the inside part is too small, a simple model results but its value to describe or predict the behaviour of the process is too restricted or too inaccurate. If the inside part is too large, a very complex model can result which is too complicated to deal with. Both situations are not attractive and have to be avoided. So, the model should contain only those components or functions that are really important for our study. A consequence will be that there are many different models for the same process, each model for its own application. Moreover, any model will always be “smaller” and describe less details than the real-life process, so
Model ⊂ Process Another important factor deals with the time range. If we want to study physiological phenomena the time scale of the processes can be very fast (milliseconds for electrical activities) to a billion years (dealing with evolution), so from 10-3 to 1016 seconds, a range of 1019 between smallest and largest time interval. This is too large! Reasonable models have a range of maximal 104. A larger range will introduce complexity, large calculation times and not comparable results. In general, processes that are faster than the selected time range of the model are assumed to have reached an equilibrium, so are assumed to be constant. Process that are slower that the selected time range of the model are so slow that they are considered too as constant in the time considered by the model. Example In studying the changes of the lung during evolution, the time range has to be at least 100.000 years or more. Then, describing in the same model also the periodic breathing of 10 seconds is both unrealistic and technically almost impossible.
Fast model
10-6
10-3
Slow model
1
103
106
109
1012
1015
1018
time
[sec]
In general we can distinguish between two different kinds of models, namely • structural or mechanistic models that have a close resemblance with the underlying process, for example, electric circuits, drawings of hydraulic or mechanical systems. Real physical components and their interconnections can be distinguished. The internal variables are clearly visible. These models incorporate all available knowledge. They are quite complex with many, sometimes unknown, parameters. • functional models that just reflect the external behaviour of a process among its inputs and outputs. The internal construction or operation is of less importance. Main goal is to represent its external behaviour. The internal processes can be hidden. Examples of functional models are mathematical models such as differential equations, state space models, transfer functions, impulse responses or their graphical representations as block diagram. Functional models are more abstract, and, consequently, more generally applicable. Mathematical models can be described by equations. For our purpose we distinguish between algebraic equations, differential equations and state space equations: In algebraic equations there is an instantaneous relation among the variables, for example • the relation between force f and displacement x of a spring: f = k.x •
the relation between pressure difference ∆p and flow φ in a pipe:
3
φ = c. ∆p
In differential equations the time behaviour or dynamics of variables is important. For example, when applying an input to a dynamic process, it takes some time before the effect becomes visible and reaches its maximum effect. Sometimes it disappears again. Such processes are called dynamic. Examples: drinking one or many glasses of beer results in a gradually increase of alcohol in the blood and after some hours in a “hang-over”. Many hours later all effects will disappear. A charged capacity can be discharged by a resistor, but it takes some time before voltage has dropped below a certain value. Pressing the gas pedal in a car will result in an increase of velocity, but not immediately the final speed will be reached. Differential equations are appropriate tools to describe dynamical phenomena. An example of a differential equation is
&y&&(t ) + 2 &y&(t ) + 5 y& (t ) + 10 y (t ) = 4u& (t ) + 8u (t ) In general we call the input variable u(t) and the output variable y(t). State space equations are also used to describe the time behaviour or dynamics of processes. Differential equations can be rewritten as state space equations and vice versa. For example the previous differential equation yields the following state space model. Besides the input u(t) and the output y(t), the internal variable x(t), the state, is introduced. An n-th order differential equation will yield a state space model with an n-dimensional state vector x(t) and n first-order differential equations, one for each of the n elements of x(t). Moreover, there is an algebraic output equation to define the output y(t) as a function of the state x(t) and input u(t).
x& (t ) = f ( x(t ), u (t )) y (t ) = g ( x(t ), u (t )) For linear functions f and g the state equations can be formulated with vectors and matrices, for example:
⎛ x&1 (t ) ⎞ ⎛ − 2 − 5 − 10 ⎞ ⎛ x1 (t ) ⎞ ⎛ 1 ⎞ ⎟ ⎜ ⎟ ⎟⎜ ⎟ ⎜ ⎜ 0 0 ⎟.⎜ x 2 (t ) ⎟ + ⎜ 0 ⎟.u (t ) ⎜ x& 2 (t ) ⎟ = ⎜ 1 ⎜ x& (t ) ⎟ ⎜ 0 1 0 ⎟⎠ ⎜⎝ x3 (t ) ⎟⎠ ⎜⎝ 0 ⎟⎠ ⎝ 3 ⎠ ⎝ ⎛ x1 (t ) ⎞ ⎟ ⎜ y (t ) = (0 4 8).⎜ x 2 (t ) ⎟ + (0).u (t ) ⎜ x (t ) ⎟ ⎝ 3 ⎠ This state equation expresses the same dynamical behaviour between u(t) and y(t) as expressed by the differential equation 1.1.1 Differential equations In differential equations not only the value of a variable y(t) is important but also its rate of change, expressed as:
dy (t ) dt
= y& (t ) = lim
∆t → 0
y (t + ∆t ) − y (t ) ∆t
This rate of change of y(t) is the time derivative of y(t), so the slope of y(t) when it is drawn as a function of time t. Besides this first-order derivative y& (t ) also higher-order derivatives of y(t) can occur in a differential equation, such as the second-order &y&(t ) , third-order derivative &y&&(t ) or n-th order derivative y(n)(t). The order n of a differential equation is determined by the highest derivative of y(t) in the equation. For calculating a unique solution of an n-the order differential equation, the values of n initial conditions have to be known, so y(0), y& (0), &y&(0),.. , y(n-1)(0). Differential equations turn out to be valuable mathematical descriptions for studying, analysing and synthesizing dynamical processes.
4
Examples: First-order differential equations with 1 initial condition and their solution Differential equation Initial condition Solution y(t) t≥0 y(0) = y0 y(t) = y0 y& (t ) = 0
y& (t ) = a y& (t ) + ay(t) = 0 y& (t ) + ay(t) = b
y(0) = y0
y(t) = y0 + at
y(0) = y0
y(t) = y0. exp(-at);
y(0) = y0
y(t) = (b/a) {1-exp(-at)} + y0.exp(-at);
Second-order differential equations with 2 initial condition and its solution &y&(t ) + 3 y& (t ) + 2 y (t ) = 1 y(0) = y& (0) =0 y(t)=1 - 2exp(-t) + exp(-2t) 1.1.2 State space equations In dynamic models we can define a state vector, besides the input and output. This state is an internal variable of the model, and represents, in general, the physical buffers of a process. If we know the state x(t) at time t=0, so x(0), and the value of the input u(t) for t≥0, then both x(t) and y(t) for t≥0 are fully defined by the state (space) equation:
x& (t ) = f ( x(t ), u (t )) y (t ) = g ( x(t ), u (t )) These state equation define the first time derivative x& (t ) of the state x(t) as an algebraic function of ONLY the state x(t) itself and the input u(t). The output y(t) is defined as an algebraic function of ONLY the state x(t) and the input u(t). When the model is of n-th order, the state x(t) will be a vector with n elements. When the model has one input and one output, both u(t) and y(t) are scalar functions. For obtaining a unique solution x(t), and so y(t), there has to be an initial condition x(0) for t = 0. Vector x(0) has n elements. When both functions f(x,u) and g(x,u) are linear in both x and u, linear state equations result:
x& (t ) = Ax(t ) + Bu (t ) y (t ) = Cx (t ) + Du (t ) Vectors and matrices: State space models utilise vectors and matrices. These are nice an attractive components from linear algebra to simplify variables and equations. Vectors have n elements, either ordered as a row (x1, x2, x3, .. , xn) or as a column (x1, x2, x3, .. , xn)T. A matrix has, for example, n rows and m columns. The dimension of a row vector is (1*n), a column vector has dimension (n*1) and this matrix has dimension (n*m). Example. The matrix equation y =Ax represents the following equations when y is a 2-dimensional and x a 3dimensional column vector:
⎛ y1 ⎞ ⎛ a11 ⎜⎜ ⎟⎟ = ⎜⎜ ⎝ y 2 ⎠ ⎝ a 21
a12 a 22
⎛ x1 ⎞ a13 ⎞⎜ ⎟ ⎟⎜ x 2 ⎟ a 23 ⎟⎠⎜ ⎟ ⎝ x3 ⎠
The associated equations are:
y1 = a11 x1 + a12 x 2 + a13 x3 y 2 = a 21 x1 + a 22 x 2 + a 23 x3 Clearly, y = Ax is a much shorter notation. Some more matrix characteristics: A.A-1 = I, with I the unity matrix. When A = BT (T from transpose), the following is valid: akl = blk.
1.2 Deriving a model 1.2.1 SI units 5
When dealing with physical systems it is a necessity to take care of appropriate units for all variables and parameters of components. Much confusion and even errors can be avoided if proper usage is made of these units. In an equation the combined effect of all units has to be the same for all individual parts that are added or subtracted. For example, in the following equations the third term is not correct, just looking at the proper units of the individual terms:
u3 = R.i + u3 + L.i R3C1 = R2 C 2 +
C3 R2 R1 R1 + C2 R3 R2
There is an international scientific standard of units (SI: Systeme Internationale). Unhappily, the USA utilizes other units (miles in stead of meters, pounds in stead of kg, etc.) which might give considerable confusion. The seven basic SI units are: physical variable length mass time electrical current temperature quantity matter light flow
symbol l m t i T n I
[m] [kg] [s] [A] [K] [mol] [cd]
unit meter kilogram second Ampere Kelvin mol candela
1.2.2 Analogies among different physical fields In several physical fields the same type of variables and components can be recognized. These variables always occur as pairs. They are selected such that their product yields power [W]. The interaction among physical components can best be understood by considering the exchange of power among them. Understanding one field makes understanding the other fields more easy. Table 2 describes the analogies among the fields electrical, mechanical, hydraulic, acoustic and thermal systems. The first variable is called effort or across variable (voltage, force, torque, pressure, temperature). The second one is the flow or through variable (current, velocity or flow). The product of effort and flow yields the physical quantity power P [W]. In these five physical fields or application areas the following generic components can be distinguished with their unique relation between their effort (e) and flow (f) variables: Resistance R
e=f·R
Capacity C
e& = C1 ⋅ f or e(t ) = e(0) + C1 ∫
Inductance L
f& = 1L ⋅ e or f (t ) = f (0) +
t
τ =0
f (τ )dτ
t 1 L τ =0
∫
e(τ )dτ
For electrical systems the effort is the voltage u [V] and the flow variable is the current i [A]. Then these 3 generic equations become: Electrical resistance R [Ω]
u=i·R
Electrical capacity C [F]
u& = C1 ⋅ i or u (t ) = u (0) + C1 ∫ i (τ )dτ
Electrical inductance L [H]
i& = L1 ⋅ u or i (t ) = i (0) +
t
τ =0
t 1 L τ =0
∫
u (τ )dτ
Only a resistance (R or b) dissipates power and so heats up. This thermal power can, in general, not be recovered and will leave the system. Capacity and inductance (C and L) are buffers, so they can store power [W] as energy [J] = [Ws]. The amount of energy stored in elucidated in table 2. The energy E(t) [J] stored in a C [F}: t
t
t 1 1 E (t ) = E (0) + ∫ u (τ )i (τ )dτ = E (0) + ∫ u (τ ).Cu& (τ )dτ = E (0) + C ∫ u (τ )du = E (0) + C[u 2 (t ) − u 2 (0)] = Cu 2 (t ) τ =0 2 2 0 0
6
In the same way, the energy stored in an inductance L [H] yields: t
t
1 1 E (t ) = E (0) + ∫ u (τ )i (τ )dτ = E (0) + ∫ Li&(τ ).i (τ )dτ = E (0) + L ∫ i (τ )di = E (0) + L[i 2 (t ) − i 2 (0)] = Li 2 (t ) τ =0 2 2 0 0 t
Power enters and leaves a buffer without a loss. When making a model, each buffer introduces a state. Because, if we know content of the buffers at time t=0, so x(0), and the value of the input u(t) for t≥0, then both x(t) and y(t) for t≥0 are fully defined by the state (space) equation x& (t ) = f ( x(t ), u (t ))
The number of states equals the number of buffers In table 2 additional variables with their units are introduced and defined for specific application areas such electrical, mechanical, hydraulic and thermal variables. For each of these areas the two most important (basic) variables are elucidated in table 1. These variables are selected such that their product equals power. Only the thermal area has a different interpretation.
Table 1. Basic variables Electrical
Mechanical translation Voltage: u [V] Force: f [N] Current: i [A] Velocity: v [m/s]
Mechanical rotation
Hydraulic/ acoustic
Thermal
Torque: M or τ [Nm] Angular speed: ω [rad/s]
Pressure: p [N/m2] Volume flow: φ[m3/s]
Temperature: T [K] Heat flow: q [W]
Table 2. Analogies among different physical fields Electrical Effort Flow
u [V] i [A]
P=u·i Power P [W] Resistance R [Ω] (damper) u=i·R
Mechanical translation f [N] v [m/s]
Mechanical rotation M or τ [Nm] ω [rad/s]
Hydraulic/ acoustic
Thermal
p [N/m2] φ[m3/s]
T [K] q [W]
P=f·v
P=τ·ω
P=p·φ
P = q !!
b [Ns/m] f=v·b
b [Nms/rad] τ=ω·b
RH [Ns/m5] lin: p = φ· RH non-lin: p = RH’ · φ2
RT [K/W] T = q · RT
φ =k p
R Capacity (spring)
C
C [F]
u& = ⋅ i 1 C
(i = iin – iout )
k [N/m] f& = k ⋅ v (f=k·x)
CH [m5/N]
k [Nm/rad] τ& = k ⋅ ω (τ=k·θ)
p& =
CT [J/k]
⋅φ
1 CH
T& = C1T ⋅ q
(φ = φin – φout )
~~
Inductance L [H] &i = 1 ⋅ u (mass, L inertia, inertance)
m [kg] v& = m1 ⋅ f
J [kgm2]
(f=m·a)
(τ = J · α )
m → v, x
J → ω, θ θ& = ω
( q = qin – qout )
LH [kg/m4]
ω& = 1J ⋅τ
φ& =
1 LH
⋅p
L Buffer with C → u state u&C = C1 i variable L→i
x& = v
v& = &x& =
i&L = L1 ∆u Energy stored in C- ½C.u2 buffer [J] Energy stored in L- ½L.i2 buffer [J]
1 m
∑ fi
ω& = θ&& =
CH → p 1 J
∑τ i
p& =
1 CH
T& = C1T
LH → φ
φ& =
1 LH
½k.x2
½k.θ 2
½CH.p2
½m.v2
½J.ω2
½LH.φ2
7
∑φ
CT → T
⋅ ∆p CT.T
∑q
i
Besides the components R, C and L there are more complex components that take care of power conversion within or between two fields, for example transformers, motors and pumps. When these components are ideal, there is no power loss. Consequently, the power at one side has to be equal to the power at the other side. These components are illustrated in table 3. Table 3. Power conversion On both sides the same variables (effort or flow) E/E transformer u2 = α · u1, i1 = α · i2 Mr/Mr gear box τ2 = α · τ1, ω1 = α · ω2 Mt/H pump f = A ( p2 – p1 ), φ = A · v On both sides different variables: E/Mr DC-motor τ = α · i, u = α · ω
P1 = u1 · i1 = u2 · i2 = P2 P1 = τ1 · ω1 = τ2 · ω2 = P2 PM = f · v = ( p2 – p1 ) · φ = PH PE = u.i = αω.τ/α = ω.τ = PM
α
E/Mt linear motor/ loudspeaker Mr/H centrifugal pump
f = α · i, u = α · v
PE = u.i = αv.f/α = v.f = PM
p2 – p1 = α · ω , τ = α · φ
PH = ( p2 – p1 )φ = αω·τ / α = ω·τ = PM
α
1.2.3 Deriving a model in 4 steps
1a: Define input u 1b: Define state x 2: Define state equation
Inputs are external variables that cannot be influenced Each buffer is a state Each state variable has its own state equation (table 1)
3: Remove all variables in right hand side, except for u and x
Use additional equations for this removal. Verify signs of state variables in equations
4: Define output equation, only as function of u and x
Outputs are given by the problem formulation or application
x& (t ) = f ' ( x(t ), u (t ), z (t ),..) x& (t ) = f ( x(t ), u (t )) y(t ) = g ( x(t ), u (t ))
A state space model can always be found by using the following 4 steps: 1. Define inputs u and determine the buffers. Each buffer has a state variable xi associated with it. Define this state. The number of independent states becomes the order of the state space equation. 2. Each state variable has its own state equation (see table 2). 3. Remove/eliminate all variables that are no inputs u and no states x from the right hand side. This goal can always be achieved by using additional equations, for example applying Kirchhoff’s laws, component equations, etc. Verify signs of state variables in right hand side. 4. Define the output equation as function of only the input and the state. 1.2.4 Sign rule The state variables at the right hand side of the state equation that correspond with the state variable at the left side will have ALWAYS a – sign. This rule is a characteristic of all physical (passive) processes. Remark: When an input u is used that depends on the state variables, this rule does not apply. The sign of the state variable in the input is then uncertain. Example: A mechanical mass-spring system yields: b
k
m&x& + bx& + kx = f or 1 1 v& = ( f − bx& − kx) = ( f − bv − kx) . m m
m f
8
x, v
1.3 Models in different physical fields 1.3.1 Electrical processes * assign a voltage (+/-) and current (direction) to each component * Each capacity Ci yields a state ui * Each inductance Li yields a state ii * Each junction of “wires”: apply Kirchhoff’s curent law: n
net current in any connection is zero or
∑i j =1
j
=0
* Each loop of components: apply Kirchhoff’s voltage laws n
net voltage along a loop in a circuit is zero or
∑u j =1
j
=0
Example A parallel connection of a resistor R [Ω], an inductance L [H] and a capacitor C [F] with a current source I [A], yields a model with input I and states uc [V] and current iL [A]: I
u 1 ( I − iL − C ) C R 1 i&L = u C L u& C =
C
R
L
If there is a series connection of a voltage source U [V] with a resistor R [Ω], an inductance L [H] and a capacitor C [F] the associated model becomes with input U and states uc [V] and current iL [A]:
1 iL C 1 i&L = (U − i L R − u C ) L u& C =
L U
R C
1.3.2 Mechanical processes translation: * Assign a velocity v (direction) and position x (direction) to each mass. * Each mass yields a state for both its velocity v and its position x. * Search for all forces in one direction for each mass. Take care of the signs. The net sum of these forces yields the acceleration in that direction. If there are more directions (x, y in 2D and x, y and z in 3D) repeat procedure for next direction.
Sometimes, the masses and accelerations are neglected. Only velocity and position are considered as relevant. Then the state equations reduce to
x& =
1 1 ( f − kx) , or, without spring constant, x& = f b b
rotation: Assign an angular velocity ω and an angle θ to each inertia. Each inertia yields a state for both its angular velocity ω and its angle θ.
Search for all torques/moments in one direction for each inertia. Take care of the signs/rotation direction. The net sum of all these torques yields the acceleration in that direction. If there are more directions (x, y in 2D and x, y and z in 3D) repeat procedure for next direction. 1.3.3 Hydraulic processes * Assign a pressure across (+/-) and flow (direction) to each component * Each hydraulic capacity CHi yields a state pi 9
* Each hydraulic inertance LHi yields a state φi The nonlinear hydraulic resistance RH, modelled with
φ = k p , is only valid for positive values of φ. If the
flow goes in the opposite direction, this resistance has to be modelled as:
φ = k .sign( p) p
.
The hydraulic capacity CH [m5/N] is defined as the ratio of the change of volume V at an increase of the pressure p:
CH =
∆V [m5/N] ∆p
For an open vessel with cross section A [m2], gravity g [m/s2] and fluid with specific mass ρ [kg/m3] and height of fluid h [m]: ∆V = A.∆h and ∆p=ρg∆h, so
CH =
A.∆h A. = . ρg∆h ρg
In stead of the pressure p, sometimes the height h is utilized as state variable:
p& =
1 CH
∑φ =
ρg
∑ φ , with p=ρgh this yields A
ρg 1 p& = ρgh& = φ , so h& = ∑ φ , ∑ A A
If the wall of a vessel or tube is elastic, that vessel/tube exhibits a compliance CH. The hydraulic inertance LH [kg/m4] represents the mass of the fluid to be accelerated by the pressure difference. The hydraulic inertance of a fluid in a pipe with length l [m] and cross section A [m2] becomes:
LH = ρ
l [kg/m4]. A
Example A vessel with cross section A [m2] contains a liquid with height h [m]. There is a flow in (φi) and a flow out (φo) of the vessel. The out flow leaves the vessel at the bottom by a nonlinear restriction k. Input: fi, state p. The model becomes:
p& =
φi
1 1 ρg (φ in − φ out ) = (φ in − k p ) = (φ in − k p ) . CH CH A
A p
When the height is selected as state, the model becomes:
1 1 h& = (φ in − φ out ) = (φ in − k ρgh ) . A A
h
k φ0
1.3.4 Thermal processes Each thermal capacity CT yields a state Ti Friction and a resistance introduce losses which are the input q [W] for a thermal process. This thermal power cannot be transformed without losses to mechanical or electrical power. Example Suppose a baby with a metabolism of qb [W] in an incubator with thermal capacity CT and resistance to environment RT and heating qh [W]. Temperature environment equals Te [K]. With inputs qh, qb and Te, state T of the incubator, the model becomes:
T − Te 1 1 (qin − q out ) = (qb + q h − ). T& = CT CT RT
10
CT
TE
qb
T RT
qh
1.3.5 Mass flow, compartment models (pharmacy kinetics)
Relevant quantities are mass m and mass flow f to be expressed in [mol] or [kg] and [mol/s] or [kg/s]. Mass is the state variable, so a buffer that receives and releases mass. Basic law: Preservation of mass in a compartment:
m& = f in − f out .
When there are many compartments with mass flows among them, the following variables are to be distinguished: mi : mass in compartment i [mol] or [kg] fij: mass flow from j to i [mol/s] or [kg/s] foi: mass loss/excretion from compartment i [mol/s] or [kg/s] ui: external supply in compartment i [mol/s] or [kg/s] Mass exchange between two compartments and environment:
m& 1 = f in − f out = f12 + u1 − f 21 − f 01 . m& 2 = f in − f out = f 21 + u 2 − f 12 − f 02 The flows fij can follow several laws: Linear fij = kij*mj
Concentrations fij = cij*yj = (kij*Vj)*(mj/Vj)
Langmuir fij = kij*(1-σi.mi)*mj
Michaelis-Menten fij = a*mj/(b+mj)
Two compartment model with linear flows:
m& 1 = −(k 21 + k 01 )m1 + k12 m2 + u1 m& 2 = k 21m1 − (k12 + k 02 )m2 + u 2
or
.
k12 ⎞ ⎛ m1 ⎞ ⎛ 1 0 ⎞ ⎛ u1 ⎞ ⎛ m1 ⎞ ⎛ − (k 21 + k 01 ) ⎟.⎜ ⎟ + ⎜ ⎜⎜ ⎟⎟ = ⎜⎜ ⎟.⎜ ⎟ . − (k12 + k 01 ) ⎟⎠ ⎜⎝ m2 ⎟⎠ ⎜⎝ 0 1 ⎟⎠ ⎜⎝ u 2 ⎟⎠ k 21 ⎝ m2 ⎠ ⎝ Example: Suppose a vessel with a constant volume V [m3] and flow fi [m3/s] in and fo = fi out. The concentration in the input flow of some material is ci [kg/m3]. Then, the mass flow into the vessel is fi*ci [kg/s]. The state is the mass m [kg] in the vessel. The concentration cv in the vessel becomes cv=m/V. Consequently, the concentration of the output flow co = cv. Inputs: fi, ci and fo , state: m. Now the state equation becomes:
m& = f i ci − f o co = f i ci − f o
m . V
m
V
fi, ci
fo
Example:
Suppose a vessel with a changing volume V = A.h [m3] and flow fi [m3/s] in and f o = k h . The concentration in the input flow of some material is ci [kg/m3]. Then, the mass flow into the vessel is fi*ci [kg/s]. The states are the mass m [kg] in the vessel and the height h [m] of the liquid. The concentration cv in the vessel becomes cv=m/V=m/(Ah). Consequently, the concentration of the output flow co = cv. Inputs: fi and ci, states: m and h. Now the state equations become: fi, ci
m& = f i ci − f o co = f i ci − f o
m A.h
A m
1 1 h& = ( f i − f o ) = ( f i − k h ) A A
h k fo
11
Example: Two cells are connected. There is an exchange of mass between both cells fij = k*(cj - ci) with ci = mi/Vi the concentration in cell i, and there are additional inputs ui and outputs f0i [kg/s] for each cell. Inputs: u1, u2, f01 and f02, states: m1 and m2. Then the state equations for the states m1 and m2 become
m& 1 = u1 − k * (
f2, c2
m1 m2 ) − f 01 − V1 V2
m& 2 = u 2 + k * (
u1 fo1
m1 m2 ) − f 02 − V1 V2
V 1 m1
V2 m2
u2 fo2
f1, c1
In this example the exchange of mass between the two volumes/buffers is driven by the difference in concentration. The larger this difference, the large the mass flow will be. The additional inputs ui and outputs f0i are modelled as mass flows, independent of concentrations. If, in a new problem, the outflows are defined as volume flows f [m3/s], then the modified model becomes
m& 1 = u1 − k * (
m1 m2 m − ) − f 01 1 V1 V2 V1
m& 2 = u 2 + k * (
m1 m2 m − ) − f 02 2 V1 V2 V2
1.3.6 Chemical reactions
State variable: amount of material X: [X] is number of mol’s of material X: [mol]
A+ B
k1 ⎯⎯→
←⎯⎯ k2
C
[A], [B], [C]: amount of mol [mol] k1 reaction constant [mol-1s-1] k2 reaction constant [s-1]
[ A& ] = [ B& ] = −[C& ] = −k1 [ A] * [ B] + k 2 [C ] Equilibrium:
[ A] * [ B ] k 2 [ A& ] = [ B& ] = [C& ] = 0 , so k1 [ A] * [ B ] = k 2 [C ] or = = constant. [C ] k1
Suppose at t = 0 [A] = a0 [mol], [B] = b0 [mol] and [C] = c0 [mol]. Define x [mol] as the increase of C, so also the decrease of A and B. Then is valid: [A](t) = a(t) = a0-x(t), [B](t) = b(t) = b0-x(t) and [C](t) = c(t) = c0+x(t). There is only one state variable because the 3 materials are tied together by the chemical reaction equation. State equation of this chemical reaction:
x& = k1 (a0 − x)(b0 − x) − k 2 (c0 + x) Because there can be no negative amount of materials, there are bounds for x: x≤a0, x≤b0 and x≥-c0. Consequently –c0≤x≤min(a0, b0). The equilibrium for x is found when x& = 0 . So, the equilibrium xe for x, and so for A, B and C, is determined by the solution of
k1 (a0 − xe )(b0 − xe ) − k 2 (c0 + xe ) = 0. Consequently, ae = a0 – xe, be = b0 – xe and ce =
c0 + xe. For the more complex chemical reaction
2A + B
k1 ⎯⎯→
←⎯⎯ k2
C+D
follows: 12
½[ A& ] = [ B& ] = −[C& ] = −[ D& ] = −k1[ A] * [ A] * [ B] + k 2 [C ] * [ D] Suppose at t = 0 [A] = a0 [mol], [B] = b0 [mol], [C] = c0 [mol] and [D] = d0 [mol]. Define x [mol] as the increase of C, so also the decrease of A and B and increase of D. Then is valid: [A](t) = a(t) = a0-2x(t), [B](t) = b(t) = b0-x(t), and [C](t) = c(t) = c0+x(t) ), and [D](t) = d(t) = d0+x(t). There is only one state variable because the 4 materials are tied together by the chemical reaction equation. State equation of this chemical reaction:
x& = k1 (a0 − 2 x) 2 (b0 − x) − k 2 (c0 + x)(d 0 + x) Because there can be no negative amount of materials, there are bounds for x: x≤½ a0, x≤b0, x≥-c0 and x≥d0. Consequently –min(c0, d0) ≤ x ≤ min(½ a0, b0). The equilibrium for x is found when x& = 0 . So, the equilibrium xe for x, and so for A, B, C and D, is
determined by the solution of k1 ( a0 − 2 x) (b0 − x ) − k 2 (c0 + x)(d 0 + x ) = 0 . Consequently, ae = a0 – 2xe, 2
be = b0 – xe, ce = c0 + xe and de = d0 + xe Example A blood vessel transports nutrient A along an organ. The blood flow is f [m3/s] and the concentration of A c1 [mol/m3]. When bloodvessel passes the organ there is an exchange of the nutrient by diffusion with diffusion constant kd [m3/s]. The volume of the blood vessel close to the organ is Vb [m3] and of the organ the volume equals Vo. In the organ a chemical reaction takes A with B to produce C with reaction constants k1 [1/(mol.s)] and k2 [1/s]. Assume [B] [mol] and [C] [mol] to be constant. Inputs f, c1, [B] and [C], states mb [mol] and mo [mol], the number of mols in the volumes Vb and Vo respectively. The model becomes
m& b = f * c1 − k d * ( m& o = k d * (
mb m o m − )− f b Vb Vo Vb
mb m o − ) − k1 [ A] * [ B ] + k 2 [C ] Vb Vo
1.4 Linearization of non-linear state equations For a nonlinear algebraic function f(x) a linear approximation δf(δx) can be calculated in some selected operating point xo. It will be the tangent df(x)/dx in that operating point {xo, f(x0)}. So, the linear approximation δf(δx) in {xo, f(x0)} becomes with x = x0 + δx:
f ( x ) = f ( x 0 + δx ) ≈ f ( x 0 ) +
δf (δx) =
∂f ( x) ∂x
∂f ( x) ∂x
( x − x0 ) = f ( x0 ) + δf (δx) x = x0
δx x = x0
This approximation is accurate when
x ≈ x0 so for small values of δx.
13
8 f(x)
7 6 5 4
f(x0)
3 2 1 0
x0
-1 -0.5
0
0.5
1
x
1.5
2
Figure 2. Nonlinear function and its linear approximation on x0. Example If f(x)= 2x2+5x+7, then the linear model for x=4 becomes
δf (δx) =
∂f ( x) ∂x
x=4
δx = (4 x + 5) x = 4 δx = 21δx
This linearisation of an algebraic function can also be applied for the linearisation of nonlinear state equations. Applying the following 3 steps yields a linear state space model 1. Define operation point (x0, u0) for non-linear system
x& = f ( x, u ) y = g ( x, u )
by assuming rest/steady state. So, x& = 0 • • Select the rest/steady state/average value of u = u0. If u(t)=100+3sin(5t), u0=100. Calculate the operating point x=x0 from 0 = f (x0 ,u0 ) . Calculate the value of y in the operating point: y0=g(x0,u0). 2. Linearise the non-linear functions f(x,u) and g(x,u) using the Taylor approximation. Assume that u=u0+δu, x=x0+δx and y=y0+δy, so the model is valid in the neighbourhood of (u0, x0 and y0) with small variations of δu, δx and δy.
x& = f ( x, u ) = f ( x0 + δx, u0 + δu ) = f ( x0 , u0 ) + .
Observe x& = ( x0 + δx ) = δx& ,
∂f ( x0 , u ) ∂f ( x, u0 ) ( x − x0 ) + (u − u0 ) ∂u u = u 0 ∂x x = x0
f ( x0 , u 0 ) = 0 , δx= x - x0 and δu= u - u0. Then we have the following
linear model in the neighbourhood of (u0, x0):
δx& =
∂f ( x0 , u ) ∂f ( x, u0 ) δx + δu ⇒ δx& = Aδx + Bδu ∂u u = u 0 ∂x u = u 0
3. In the same way the linearised output equation becomes:
14
y = y0 + δy = g ( x, u ) = g ( x0 + δx, u0 + δu ) = g ( x0 , u0 ) + Observe g ( x0 , u 0 )
∂g ( x, u0 ) ∂g ( x0 , u ) ( x − x0 ) + (u − u0 ) ∂x ∂ u x = x0 u =u0
= y 0 , δy= g(x,u) - y0, δx= x - x0 and δu= u - u0, then we have the following linear
model in the neighbourhood of (u0, x0):
δy =
∂g ( x, u0 ) ∂g ( x0 , u ) δx + δu ⇒ δy = Cδx + Dδu ∂x ∂u u = u 0 x = x0
Remark: partial derivatives If the model is second order, then both x(t) = {x1(t), x2(t)}T and f(x,u) = {f1(t), f2(t)}T are vectors with 2 elements. Suppose there is only one input u. The partial derivatives of f(x,u) and g(x,u) with respect to x and u become:
⎛ ∂f 1 ⎜ ∂f ( x, u ) ⎜ ∂x1 = ⎜ ∂f 2 ∂x ⎜ ∂x ⎝ 1
∂f1 ∂x 2 ∂f 2 ∂x 2
⎞ ⎛ ∂f1 ⎞ ⎟ ⎜ ⎟ ∂g ( x, u ) ⎛ ∂g f x u ∂ ( , ) ⎟, = ⎜ ∂u ⎟ , = ⎜⎜ ⎟ ∂u ∂x ⎜ ∂f 2 ⎟ ⎝ ∂x1 ⎜ ⎟ ⎟ ⎝ ∂u ⎠ ⎠
∂g ∂x 2
⎞ ∂g ( x, u ) ⎛ ∂g ⎞ ⎟⎟ , =⎜ ⎟ ∂u ⎝ ∂u ⎠ ⎠
Example:
x& = x 2 − 4 .
Operating points: x − 4 = 0 , x01= 2 and x02 = -2. Derivative of f(x) with respect to x becomes 2x0 and the derivative with respect to u is zero. Consequently, the following 2 linear models result: x01 = 2 : δx& = 2 x0δx + 0δu = 4δx , which is an unstable linear model 2
x02 = -2 :
δx& = 2 x0δx + 0δu = −4δx , which is a stable linear model
Stability can be shown by looking at the equations. Assume for x01 δx is positive. So will its derivative, which implies that δx will become even larger and bigger. δx will not return to zero and so x not to x01. This behaviour is called unstable. Assume for x02 δx is positive. Then its derivative is negative, which implies that δx will become smaller. δx will to zero and so x to x02. This behaviour is called stable. A better explanation is given in chapter 2. Example: Second-order system with state x = (z, y)T; y is now not the output but a state variable!
z& = −2 z + zy − 10 y& = −3 y + y 2
Operating points:
⎛ ∂f1 ⎜ ∂f ( x, u ) ⎜ ∂z = ⎜ ∂f 2 ∂x ⎜ ∂z ⎝
.
z& = y& = 0 , so x0 = (z0, y0)T = (-5, 0)T or (10, 3)T. ∂f 1 ⎞ ⎛ ∂f 1 ⎞ ⎟ ⎜ ⎟ ⎛0⎞ 2 − + y z ⎛ ⎞ ∂ ( , ) f x u ∂y ⎟ ⎟⎟ and = ⎜⎜ = ⎜ ∂u ⎟ = ⎜⎜ ⎟⎟ ∂f 2 ⎟ ⎝ 0 − 3 + 2y⎠ ∂u ⎜ ∂f 2 ⎟ ⎝ 0 ⎠ ⎜ ⎟ ⎟ ∂y ⎠ ⎝ ∂u ⎠ ⎛ − 2 − 5⎞ ⎟δx with eigenvalues -2 and -3, so stable. − 3 ⎟⎠
For the first operating point (-5, 0)T δx& = ⎜⎜ ⎝ 0
⎛ 1 10 ⎞ ⎟δx with eigenvalues 1 and 3, so unstable. 3 ⎟⎠
For the second operating point (10, 3)T: δx& = ⎜⎜ ⎝0
15
2. Transfer functions Complex numbers The Laplace transformation and the frequency response utilize complex numbers. A complex number c can be described with • real (a) and imaginay (b) part: c = a + jb • modulus (|c|) and phase (ϕ = arg(c)) with
c = a 2 + b 2 and ϕ = arg(c) = atan(b/a), so c = |c|.exp(jϕ) j indicates the imaginary number with j = −1 j = − 1 c1 + c2 = (a1 + a2) + j(b1 + b2) c1 - c2 = (a1 - a2) + j(b1 - b2) c1 . c2 = (a1a2 - b1b2) + j(a1b2 + a2b1) = |c1 |.| c2| .exp(ϕ1 + ϕ2) c1/c2 = (|c1|/|c2|). exp(ϕ1 - ϕ2) 2
exp(c) = exp(a+jb) = exp(a).exp(jb) |exp(c)| = exp(a) arg(exp(c)) = b Remark: Mathemathicians (and Matlab) utilize the i to indicate the complex number. Here we use the j, to make a distinction between current i and complex number j.
2.1 Laplace Transformation The Laplace transformation assigns for a real-valued variable y(t) as function of time t a complex variable Y(s) as function of the complex variable s=σ+iω ∞
L{ y (t )} = Y ( s ) = ∫ y (t )e − st dt 0
For some functions y(t), the associated value of Y(s) is shown in the table. If Y(s) has been calculated, the associated value of y(t) can also be retrieved from this table. y(t), t≥o impulse: δ(t) step: U(t) ramp: t exp(-at) t.exp(-at)
Y(s) 1 1/s 1/s2 1/(s+a) 1/(s+a)2
Y(s) ω/(s2+ω2) s/(s2+ω2) ω/((s+a)2+ω2) (s+a)/((s+a)2+ω2)
y(t), t≥o sin(ωt) cos(ωt) exp(-at)sin(ωt) exp(-at)cos(ωt)
When both input u(t) and y(t) have been Laplace transformed into U(s) and Y(s), the relation been both becomes the transfer function H(s), with Y(s)=H(s).U(s). It turns out that it is much easier to analyse H(s) than the differential equation of y(t) and u(t). Very attractive properties of the Laplace transformation are the following relations:
L{ y& (t )} = sY ( s) − y (0) and L{&y&(t )} = s 2Y ( s) − sy(0) − y& (0) , etc. Applying these relations to a differential equation,
&y&&(t ) + 2 &y&(t ) + 5 y& (t ) + 10 y (t ) = 4u& (t ) + 8u (t ) yields
s 3Y ( s) − s 2 y(0) − sy& (0) − &y&(0) + 2s 2Y ( s) − 2sy(o) − 2 y& (0) + 5sY ( s) − 5 y(0) + 10Y ( s) = 4sU ( s) + 8U ( s)
16
4s + 8 s 2 y (0) + s( y& (0) + 2 y (0)) + &y&(0) + 2 y& 0) + 5 y (0) + U ( s ) s 3 + 2 s 2 + 5s + 10 s 3 + 2s 2 + 5s + 10
Y (s) =
The transfer function H(s) of the system is the relation between Y(s) and U(s) {H(s)=Y(s)/U(s)} neglecting all initial conditions. So, H(s) of this differential equation becomes:
H ( s) =
4s + 8 s + 2s 2 + 5s + 10 3
We call the roots of the denominator of H(s) its poles and the roots of the numerator its zeros. Once Y(s) is known, y(t) can be calculated by using partial fraction expansion. First the n roots si (poles) of denominator D(s) of Y(s) have to be calculated and next the n residues Ai
Y ( s) =
n A N (s) =∑ i D( s) i =1 s − si
with, if no multiple poles occur:
Ai = lim ( s − si )Y ( s ) s → si
The time response y(t) becomes, after transforming each separate root as an exponential function: n
y (t ) = ∑ Ai e sit i =1
If there are pairs of complex poles (s1,2 = -σ ± jω), so (s + σ)2 + ω2 = s2 + 2ζωns + ωn2 , these quadratic terms can be transformed as pair:
ω bs + c s +σ c − σb =b + 2 2 2 2 ω (s + σ ) 2 + ω 2 (s + σ ) + ω (s + σ ) + ω c − σb −σt y (t ) = be −σt cos(ωt ) + e sin(ωt )
Y ( s) =
ω
Examples:
y& (t ) + ay(t) = bu(t) with y(0) = yo and u(t) the unit step. Then y y y A A b b + o = 1+ 2 + o Y (s) = U (s) + o = s + a s( s + a) s + a s s+a s+a s+a A1 = s
b s( s + a)
s =0
=
b b and A2 = ( s + a ) a s( s + a)
s=− a
=
b −a
y(t) = (b/a){1-exp(-at)} + yo.exp(-at) for t ≥ 0. The influence of several parameters (σ,ω,ζ,ωn) of the transfer function H(s) concerning (complex) roots of the denominator (poles) can be linked with the parameters (t s, t p, M p, t r) of the step response of the system.
17
parameter real part imaginary part angle modulus
•
symbol -σ ω ζ=σ/ωn 0≤ζ≤1 ωn
name absolute damping frequency relative damping natural frequency
influence stability, settling time t s peak time t p overshoot M p rise time t r
Design real part si< 0, σ>0, large σ large ω large ζ, e.g. ζ>0.7 large ωn
The absolute damping σ indicates the speed with which the response settles. For positive values of σ (so negative roots!) the response will reach a constant value. For t=3/σ, exp(-σt)=exp(-3)=0.05. Then y(t) is about within 5% of its final value. For t=5/σ, exp(-σt)=exp(-5)=0.007. Consequently, σ determines the settling time t s of a response. The larger σ, the faster, the better. The imaginary part, ω, determines the frequency in y(t). The higher ω, the higher the frequency and also the smaller the peak time t p will be. The relative damping ζ {ζ=sin(θ)} determines the overshoot M p of a step response. For ζ=0 , (θ=0), there is no damping and the overshoot will be 100%. For ζ=1, (θ=90o) there are no oscillations (no imaginary poles and so no overshoot). A design value is for example ζ=0.7 which yields about 5% overshoot. The natural frequency ωn determines the rise time of a response. The higher ωn, the faster the response and the smaller the rise time t r.
• •
•
1.4
3 Im(s)
2.5
1
vss
ζ, θ
2
ω
0.8
1.5
0.6
ωn
1
0.4 0.2
0.5
0
0 -1.5
1.2
v(t) Mp
σ -1
Re(s) -0.5
-0.2
0
tp 0
t
ts 2
4
6
8
10
Figure 3. Relation between parameters of H(s) and step response. Initial value theorem Without using the inverse Laplace transformation, it is possible to calculate the value of y(t) just after t=0, indicated with y(0+), directly from Y(s):
y (0 + ) = lim y (t ) = lim sY ( s ) t →0
s →∞
Final value theorem The final value, or steady state value yss, of y(t) can be calculated for stable models:
y (∞ ) = y ss = lim y (t ) = lim sY ( s ) t →∞
Example: Suppose
Y (s) =
H ( s) =
s →0
1 4s + 8 and u(t) is a unit step, so U ( s ) = , so 2 s s + 2s + 5s + 10 3
4s + 8 , then y(0+) = 0 and y(∞) = 0.8. 2 s ( s + 2 s + 5s + 10) 3
MATLAB: input numerator polynomial:
num = [4 8]; 18
input denominator polynomial: den = [1 2 5 10]; roots of denominator: root = roots(den); eigenvalues of matrix B eigenvalues = eig(B); residues [A, p, k] = residue(num, den); A: residues, p: poles (roots), k constant if degree(num)=degree(den) Matlab calculates complex residues for complex poles. Taking two conjugated poles together, as a quadratic term, yields real coefficients and is to be preferred. hs = tf(num, den); transfer function step response yt = step (hs); draw step response plot(yt) Laplace (symbolic) laplace(exp(-2t)) yields 1/(s+2) inverse Laplace (symbolic) ilaplace(1/(s+2)) yields exp(-2t)
2.2 Relations among differential equation, state space model and transfer function 2.2.1 Derivation of a state space model from a differential equation Suppose the following differential equation is available for the description of a process:
&y&&(t ) + 2 &y&(t ) + 5 y& (t ) + 10 y (t ) = 4u& (t ) + 8u (t ) Then, the associated transfer function can be derived:
Y ( s) =
4s + 8 U (s) s + 2s 2 + 5s + 10 3
Introduce the intermediate variable Z(s) and distinguish between the numerator N(s) and the denominator D(s) of H(s). Then, we have
1 U ( s) D( s ) 1 1 Z (s) = U ( s) = 3 U ( s) 2 D( s ) s + 2 s + 5s + 10
Y ( s ) = H ( s )U ( s ) = N ( s )
Z ( s ){s 3 + 2 s 2 + 5s + 10} = U ( s ) or &z&&(t ) + 2&z&(t ) + 5 z& (t ) + 10 z (t ) = u (t ) Y ( s ) = N ( s ) Z ( s ) = (4 s + 8) Z ( s ) There are many possibilities to define the state. One way is to select x(t)=(x1(t), x2(t), x3(t))T with
x3 (t ) = z (t ) x 2 (t ) = z& (t ) = x& 3 (t ) x1 (t ) = &z&(t ) = x& 2 (t ) Then we have
x&1 (t ) = &z&&(t ) = −2&z&(t ) − 5 z& (t ) − 10 z (t ) + u (t ) = −2 x1 − 5 x 2 − 10 x3 + u (t ) x& 2 (t ) = x1 (t ) x& 3 (t ) = x 2 (t ) Moreover, the numerator uses Y(s)=N(s)Z(s)=(4s+8)Z(s) or
y (t ) = 4 z&(t ) + 8 z (t ) = 4 x2 (t ) + 8 x3 (t ) Combining both equations yields the following state space model:
19
⎛ x&1 (t ) ⎞ ⎛ − 2 − 5 − 10 ⎞ ⎛ x1 (t ) ⎞ ⎛ 1 ⎞ ⎟ ⎜ ⎟ ⎟⎜ ⎟ ⎜ ⎜ 0 0 ⎟.⎜ x 2 (t ) ⎟ + ⎜ 0 ⎟.u (t ) ⎜ x& 2 (t ) ⎟ = ⎜ 1 ⎜ x& (t ) ⎟ ⎜ 0 1 0 ⎟⎠ ⎜⎝ x3 (t ) ⎟⎠ ⎜⎝ 0 ⎟⎠ ⎝ 3 ⎠ ⎝ ⎛ x1 (t ) ⎞ ⎜ ⎟ y (t ) = (0 4 8).⎜ x 2 (t ) ⎟ + (0).u (t ) ⎜ x (t ) ⎟ ⎝ 3 ⎠ The same equations can also be selected for drawing a block scheme (figure 4), a graphical representation of the (mathematical) equivalent descriptions of differential equation, the transfer function and the state space model. Y 4
8
x1
U
Z 1/s
1/s
1/s x2
x3
-2 -5 -10 Figure 4. Block scheme of a mathematical model Remark: The coefficients of the denominator are visible (with a minus sign!) on the first row of the A matrix. The coefficients of the numerator are placed in the C-matrix. This state space model is just one of the many possible state space models from a differential equation or transfer function.
MATLAB [A, B, C, D] = tf2ss([4 8], [1 2 5 10]);
⎛1⎞ ⎛ − 2 − 5 − 10 ⎞ ⎜ ⎟ ⎜ ⎟ 0 0 ⎟, B = ⎜ 0 ⎟, C = (0 4 8), D = (0 ) A=⎜ 1 ⎜0⎟ ⎜ 0 1 0 ⎟⎠ ⎝ ⎠ ⎝ Matlab calculates a different, but also correct, state space model! 2.2.2 Derivation of the transfer function from a state space model
x& (t ) = Ax(t ) + Bu (t ) y (t ) = Cx (t ) + Du (t ) Exercise the Laplace transformation on both sides and neglect initial conditions (x(0)=0):
x& (t ) = Ax(t ) + Bu (t ) sX ( s) = AX ( s) + Bu ( s ) X ( s) = ( sI − A) −1 BU ( s) Y ( s ) = CX ( s ) + DU ( s) = {C ( sI − A) −1 B + D}U ( s) or
H ( s) = C ( sI − A) −1 B + D 20
MATLAB [num, den] = ss2tf(A, B, C, D); 2.2.3 Derivation of the differential equation from a transfer function If the transfer function H(s) is available, then Y(s)=H(s).U(s). For example:
Y ( s) =
4s + 8 U (s) s + 2s 2 + 5s + 10 3
Multiply both sides with the denominator of H(s) yields:
s 3Y ( s) + 2s 2Y ( s) + 5sY ( s) + 10Y ( s) = 4sU ( s) + 8U ( s) Recognize that snY(s) in the Laplace domain corresponds with y(n)(t), the n-the derivative of y(t), in the time domain. Consequently, the associated differential equation is
&y&&(t ) + 2 &y&(t ) + 5 y& (t ) + 10 y (t ) = 4u& (t ) + 8u (t ) 2.3 Frequency response Suppose Y(s) = H(s)U(s). Both Y, H, U and s are complex variables. The Laplace transformation allows the calculation of transient responses, such as the response on a step or an impulse. For time going to infinity, only the steady state response remains. In the complex Laplace operator s=σ+jω σ represents the transient. When dealing with frequency responses, only the steady state response is of interest. Consequently, σ disappears (=0) and the Laplace operator s reduces to jω, or = jω H ( s) ⎯s⎯ ⎯→ H(jω)
The frequency response H(jω) is also a complex number, with modulus H ( jω ) and phase arg(H(jω)). When the input signal of system H(s) is a sine wave ( u (t ) = uˆ sin(ω 1t ) ) then the output signal becomes
y (t ) = yˆ sin(ω 1t + φ ) with yˆ = H ( jω 1 ) .uˆ and ϕ = arg( H ( jω 1 )). So, the modulus H ( jω ) is the “gain” of the process as function of the frequency. Clearly, this gain is frequency dependent. When the degree of the numerator is smaller than the degree of the denominator, so there are more poles than zeros, the modulus H ( jω ) decreases to zero for increasing values of the frequency. Example: If H(s) = 4/(s+4), then H(jω) = 4/(jω+4). The modulus of H yields |H(jω)| = 4/√(ω2+16). When the input signal equals u(t) = 5sin(2t), then for ω = 2 the frequency response of H(jω): |H |= 4/√(22+16) = 0.9 and arg(H) = atan(2/4) = 0.46 [rad]. Consequently, the output becomes y(t) = 0.9*5sin(2t+0.46) = 4.5*sin(2t+0.46).
MATLAB hs=tf(num, den); bode(hs)
21
3. Regulation Many physiological processes are regulated to guarantee that the physical or chemical conditions of the body remain within acceptable boundaries, in spite of many internal and external disturbances. For example the body temperature, the blood pressure and the CO2 concentration in the blood, are kept reasonably constant, in spite of (large) varying activity of the body. This constant operating point can be realized by some kind of feedback. The values of signals are measured and are used for making corrections. The basic idea of feedback is illustrated in figure 5.
y 0. 8 f(x)
0.6
f(x) x
0.4
g(y)
g(y)
0.2
0
y
x
-0.2 -0.5
0
0.5
1
1.5
2
2.5
3
Figure 5. Static regulation.
3.1 Static feedback loop The two processes y = f(x) and x = g(y) are tied together, by means of a feedback loop. Both processes determine the operating point (x0, y0), namely x0 = g{y0) and y0 = f(x0), so x0 = g{ f(x0)} or y0 = f{ g(y0)}. Example:
x = 3− y y = x −1
The “solution”, or operating point, has to satisfy both equations. Consequently x0 = 3 - (x0 - 1), so x0 = 2 and y0 = 1. The operation points (x0, y0) = (2, 1) satisfies both functions and is the cross-point of both graphs f(x) and g(y). When x and/or y are disturbed they will return to the operating point determined by the 2 functions f(x) and g(y). When these functions change, the operating point will not be the same, but shift a bit. For example, if f(x) changes from y = x - 1 into y = x - 2, the new operating point becomes (2.5, 0.5). As long as these two graphs cross each other, so have a common point (x0, y0), there will be one or more operating points. Still, having such a point does not guarantee that the operating point is stable. It is possible that a small deviation from the operating point may result in a decrease of this deviation (so convergence to the operation point) or an increase (so divergence). The first situation is called stable and the second one unstable. For a physiological process, stable operating points are required. Example A pump (heart) takes care that there is a pressure between its input and output. Owing to internal friction, the pressure drops when the (blood) flow increases. Pipes (arteries, veins) have a resistance when there is a (blood) flow. The larger the flow, the larger the pressure drop over these pipes. Suppose P is the pressure difference over the pump and over the pipes and call F the flow. Then, the pump can be modelled as P=f(F) and the pipes as F=g(P) or P=g’(F). The cross-point of both functions yields the operation point (P0,F0) that satisfies both the characteristics of the 22
pump and pipes. Suppose: pump: P = P0 - a.F2 for the pump and P = b.F2 for the pipes. Then, the operating point of this pump with these pipes will be (P0, F0) with
F0 =
P0 b L P0 = P0 a+b a+b
Both pump and pipes determine together the value of the operating pump. There is no causality that argues that the pump determines the pressure and the pipes the flow or vice versa. Both do.
3.2 Dynamic feedback loop In the previous paragraph the functions f and g are assumed to be static functions. However, these functions can also be considered being dynamic, e.g. represented by transfer functions, for example X(s) = H1(s)Y(s) and Y(s) = H2(s)X(s). Then, X(s) = H1(s){H2(s)X(s)}, or
X (s) =
H 2 ( s) 1 and Y ( s ) = 1 − H 1 ( s) H 2 ( s) 1 − H 1 ( s) H 2 ( s)
Now, the stability is fully determined by the poles of the denominator {1- H1(s)H2(s)=0)}. If they all have negative real parts, this feedback loop is stable. Deviation from the operating point will be (partly) suppressed, and the system will return to its operating point. Example: Prey - predator relationship Rabbits and foxes in the same ecological area depend on each other. This can be formulated with the following model with r the number of rabbits and f the number of foxes. Rabbits increase (a) if not too much foxes (b) and sufficient food (c) is available. Without rabbits (e) foxes will die (d). In this model r.f represent the probability that a fox meets a rabbit, which is negative for the rabbits and positive for the foxes. If there are too many rabbits and the food becomes a limitation, their growth will be restricted. This factor is modelled with a quadratic term r2:
r& = a.r − br. f − c.r 2 f& = −df + er. f For suitable values of the coefficients a, .., f, the stable operating point (besides (0,0)) can be determined:
d e ae − cd f0 = be
r0 =
Feedback loops For technical control systems, the following control scheme is often implemented, as elucidated in figure 6. It clearly recognises an input u(t), an output y(t) and a reference r(t), which represents the required value of the output. Moreover, a disturbance d(t) is added which reflects disturbances and not-modelled parts of the physical system. The process, indicated by its transfer function H(s), has input u(t) and output y(t). A controller has been added with input e(t) = r(t) - y(t). It calculates an appropriate values u(t) for the input of the process. Nearly always the control requirements are stated such that y(t) has to follow its reference r(t) “as good as possible”. So, the goal is to keep e(t) (the error) as small as possible for different values of r(t) and d(t), and maybe, also for changing values of H(s). D(s) R(s) +
U(s)
E(s)
C(s)
+
H(s)
-
Figure 6. Unit feedback loop. 23
Y(s)
Using some calculus, we obtain: Y(s) = D(s) + H(s).U(s) U(s) = C(s).E(s) E(s) = R(s) - Y(s)
Y(s) = D(s) + H(s).C(s).[R(s) - Y(s)]
or
Y (s) =
H ( s )C ( s) 1 R( s) + D( s ) 1 + H ( s)C ( s ) 1 + H ( s )C ( s)
and
E ( s) = R( s) − Y ( s) =
1 1 R( s) − D( s) 1 + H ( s)C ( s ) 1 + H ( s)C ( s )
If we can guarantee that |H(s)C(s)|>>1, the influence of D(s) will be reduced and Y(s) ≈ R(s). The only requirement we make is that the controlled system is stable and that the product of process and controller is (very) large. The larger, the better y(t) will fit r(t). So, a well-designed controller C(s) can realise that the output of a process will follow the reference value r(t) in spite of the presence of disturbances d(t) and also in spite of changes in H(s). The requirement is that |H(s)C(s)|>>1, so neither the model H(s) of the process nor the disturbance d(t) have to be accurately known. Feedback as a mechanism for life to survive Feedback is one of the more powerful characteristics of life. In spite of large variations in physiological processes (H(s)) and large variations in the environment (D(s)), proper feedback allows a cell, an organ, an organism, individuals or even populations to maintain stable values for important physiological variables such as concentration, temperature and pressure. In these more stable conditions the cell, organ, organism or individual can survive and even prosper. Differences between physiological and technical feedback loops. D(s)
f(x) x
R(s) +
C(s)
y
H(s)
+
Y(s)
g(y) Figure 7. Comparison between physiological and technical feedback loops
In a technical feedback loop typically can be recognized: o a sensor to sense or measure the variable to be controlled o a controller with a dedicated task to minimize the error between reference and measured output o an actuator that uses the controller signal as input and influences the process, e.g. a motor or a valve. In general there is only one sensor and one actuator. In physiological feedback loops there are several processes that together determine the stable operating point of these processes. In many situations no clear reference signal can be observed. The stable operating point is established by the cross-point of two processes. Really important control loops in the human body, for example for regulating the blood pressure, have many sensors, many actuators and many “controllers” to maintain stability in many, quite different situations. Even complete failure of one or more sensors or actuators can be accepted because other control loops with different sensors and actuators will replace the function of the broken ones. So, feedback or control loops in physiology turn out to be more robust and more reliable than technical ones. These technical solutions have a higher performance for their dedicated task. Outside that task or with severe disturbances they will (completely) fail. Reliability turns out to be important for live than performance.
24
Remark: In general a process has the characteristic that |H(jω)| decreases for high frequencies. So, disturbance rejection and reference following is only guaranteed for frequencies such that |H(jω)C(jω)|>>1. When a controller has a pole in (0,0), just like an integrator, the modulus or gain will become very large when the frequency nears 0. So, perfect disturbance rejection and perfect following of the reference can be expected in steady state when using this type of controllers.
3.3 Steady state error With dynamic systems we are interested in both the transients and the behaviour in steady state. In steady state all signals and variables have reached a constant value. With the final value theorem it is possible to calculate the steady state value of variables. In this section the effect of the controller C(s) on the steady state error ess will be illustrated.
e(∞ ) = e ss = lim e(t ) = lim sE ( s ) t →∞
s →0
If r(t) and d(t) are step functions, R(s)=rss/s and D(s)=dss/s, the value of ess becomes:
1 1 R( s) − D( s) 1 + H ( s)C ( s ) 1 + H ( s)C ( s ) rss d ss 1 1 } ess = lim sE ( s ) = lim sE ( s ) = lim s{ − 1 + H ( s )C ( s ) s 1 + H ( s )C ( s ) s s →0 s →0 s →0 1 1 ess = rss − d ss 1 + H (0)C (0) 1 + H (0)C (0) E ( s) = R( s) − Y ( s) =
When C(s) and/or H(s) have a pole in the origin (0, 0) of the complex plane, the product H(0).C(0) will become very large, reducing the steady state error ess to zero. In technical control loops large values of C(0) can easily be realized with PI-controllers which have a Proportional gain K and an integrator, for example:
C (s) =
a K (s + a) = K (1 + ) s s
In physiological control loops this large gain is achieved by using steep functions with a large derivative in the operating point.
25
4. Exercises / Oefenopgaven 4.1 Refreshing mathematics
&y& + 3 y& + 2 y = 6 met y (0) = 0 en y& (0) = 0 is
1. De oplossing van de differentiaalvergelijking −t
−2t
y(t ) = 3 − 6e + 3e . Bereken y& (t ) en &y&(t ) . Laat zien dat y (t ) inderdaad de oplossing is van de differentiaalvergelijking en dan y (0) = 0 en y& (0) = 0 . 2. Bereken H zoals in onderstaand schema gedefinieerd: R
Y
H
voor de volgende blokschema’s. a. H2
+
R
+
Y
H3
H1
b. +
R
H1
Y
H2
-
H3
c. H2 +
R +
-
H1
-
H3
Y
a = 2 + 3i , b = 4 + i . Bereken c1 = a + b , c 2 = a − b , c3 = a ∗ b
3. Gegeven is
Bereken
c = e a ⋅ e b met c en ϕ = arg(c) .
4. Gegeven is:
y1 = 3 x 2 + 2 x1 + 5u1 y 2 = 4 x1 + 7 x 2 − 6u 2
⎛y ⎞ y = ⎜⎜ 1 ⎟⎟ ⎝ y2 ⎠
y = Cx + Du Bereken C en D . T Bereken det(C ) en C . 5. Gegeven is:
x&1 = −2 x1 + 3 x 2 + u x& 2 = x1 − 3x2 + x3
x& 3 = 2 x1 − 5 x3 + u Bereken A en B .
26
⎛ x1 ⎞ ⎜ ⎟ x = ⎜ x2 ⎟ ⎜x ⎟ ⎝ 3⎠
⎛x ⎞ x = ⎜⎜ 1 ⎟⎟ ⎝ x2 ⎠
⎛u ⎞ u = ⎜⎜ 1 ⎟⎟ ⎝ u2 ⎠
x& = Ax + Bu
m [kg ] en zwaartekracht g [m / s 2 ] . De vallende massa ondervindt een wrijving ter grootte b ⋅ v [N ] met v de snelheid van de massa.
6. Gegeven een massa
a. Bereken
v&(t )
Op een bepaald moment vertraagt of versnelt b. Wat is dan de snelheid v ? c. Schets v(t ) als v(0) = 0 .
m
m niet meer ( → v&(t ) = 0 ).
4.2 Electrical processes 1.
Ub L C1
C2
R
Uitgang: i1 a. Geef de ingang en de toestand(en). b. Geef het toestandsmodel en de A, B matrices.
i1
2.
Uitgang: u R2 .
R2
a. Geef de ingang en de toestand(en). b. Geef het toestandsmodel en de A, B matrices.
L
R1
Ub
C2
C1
4.3 Mechanical processes 1.
f
k12
k23 m3
m2
m1
b1
b2
b12
Uitgang: x 2 . a. Geef de ingang en de toestand(en). b. Geef het toestandsmodel.
b3
b23
x2
2.
f
Uitgang: v . c. Geef de ingang en de toestand(en). d. Geef het toestandsmodel.
x,v
m b
α mg
3.
Ingang: τ ; Uitgang: ω 2 . e. Geef de toestand(en). f. Geef het toestandsmodel en de A, B, C , D matrices.
b
τ k2
k1
I1
I2
27
4.
Ingang: vermogen P ( = τ ⋅ ω ) ; Uitgang: ω . a. Geef de toestand(en). b. Geef het toestandsmodel.
P
bω2
I 5.
a. Geef de toestand(en). b. Geef het toestandsmodel.
R
I
k f
m mg 6. Bepaal ingang, toestanden en toestandsmodel van bijgaande schakeling.
k1
⎛ x1 ⎞ ⎟⎟ . ⎝ x2 ⎠ Bereken de A, B, C , D matrices.
b1
m1
Kies als uitgang ⎜⎜
x1 b2
m1g m2 m2g
k2 x2
f
4.4 Electromechanical processes 1.
i L
τ ω
R
Ud
α
Een patiënt fietst op een hometrainer. Hij levert een koppel τ [Nm] . Aan de as is een ideale dynamo gekoppeld, die een batterij oplaadt
b
I
a. Bepaal ingang, toestanden en toestandsmodel met P=i*Ud [W] als uitgang. b. Als de uitgangsspanning gelijk is aan U 0 [V ] , wat is dan het geleverde koppel toerental
ω
τ
en het
in rust (steady state) ?
2. Bepaal een toestandsmodel van het volgende proces:
Ub
k : torsieveer bi : wrijvingen/dempers
α
R
L
k b1
b2
I2
I1
28
I i : inertia
α : ideale motor
4. k
spier m2
R
α
L
m1
k,b
m1
f m2 b
Een spier is aangesloten op een lineaire motor. Aan de elektrische klemmen is een weerstand geschakeld. De spanning u over de weerstand R wordt gemeten en is de uitgang. De spier wordt beschreven met een kracht f (t ) , die via een parallel geschakelde veer k en demper b werkt op de lineaire motor. Geef uitgang, toestand(en) en toestandsmodel.
4.5 Hydraulic processes Bepaal ingang, toestanden en toestandsmodel van dit hydraulisch proces. De vloeistof heeft een
φ1
1.
h1
ρ [kg / m 3 ] en de gravitatieversnelling 2 is g [m / s ] . φ = k p . dichtheid
A1
A2 p2
p1 k12
h2
Als k 2 = 2 en k12 = 1 , bereken dan h1 en h2 in steady state.
k2 φ0
1. In een bloedvat wordt een bloedstroom
φ [m 3 / s ]
gepompt. De compliantie (hydraulische
capaciteit) van het bloedvat bedraagt C H en de relatie tussen druk en flow in het bloedvat is
φ = k p . De uitgang is de druk over dit bloedvat. Bepaal ingang, toestand(en) en toestandsmodel. 2. Het hart pompt bloed in de slagader. Deze heeft een compliantie C H . De haarvaten hebben een weerstand beschreven door druk van
φ = k p . Het bloed heeft een inertia L H . Het hart verhoogt de
Pader tot Pslagader [ N / m 2 ] .
Bepaal ingang, toestand(en) en toestandsmodel van deze eenvoudige bloedsomloop. Veronderstel de druk in de ader constant. 4.
Voor longonderzoek blaast een patiënt een ballon op. Het mondstuk heeft een pneumatische weerstand R1 (∆p = R1φi ) . De ballon lekt met weerstand R2 naar de omgeving
(∆p = R2φ0 ) . De ballon heeft een pneumatische ‘capaciteit’ C p . Er zijn 2 situaties:
a. De patiënt heeft een constante druk
pi .
b. De patiënt geeft een constante volumestroom
φi .
Bereken in beide gevallen een toestandsmodel (ingang, toestand, vergelijking) en bepaal de druk in de ballon in steady state.
29
A1
5.
k1
R Ub
α
L
Via een spanningsbron vloeistof ( ρ [kg / m en k 2 ( φ i
A2
3
k2
m
Ap g[m/s2]
U b wordt met een lineaire motor ( α ) een pomp aangedreven. Deze verplaatst
] ) tussen 2 vaten met oppervlakte A1 en A2 . Er zijn 2 hydraulisch weerstanden k1
= k i pi ).
Beschrijf het model, als vloeistof wordt gepompt van vat 1 naar vat 2.
4.6 Thermal systems 1)
Een schaatser schaatst in een isolerend pak met thermische weerstand RT naar de omgeving met temperatuur
T0 . Zijn warmtecapaciteit is CT . Zijn metabolisme levert q [W ] . Wat is zijn
temperatuur? Als door harde wind de effectieve waarde van RT halveert, hoe groot moet dan zijn metabolisme zijn om dezelfde lichaamstemperatuur te houden? 2
R1
C1
C2 R12
T0
q1
R2
Twee organen liggen tegen elkaar aan. Beide bevinden zich in een vloeistof met temperatuur T0 . Beide hebben een interne verwarming
qi en warmteweerstand Ri naar de omgeving en
R12 naar elkaar.
q2
Geef het toestandsmodel.
4.7 Pharmakinetic models 1.
g(t)
L
B
E 1
2
3
De bloedbaan wordt gedacht te bestaan uit 3 secties: Long, Body en Extremiteiten (L, B en E). Zie de figuur hiernaast.
2. Een patiënt krijgt een injectie g (t ) in zijn arm. De stof g (t ) [ kg ] wordt via de bloedbaan verspreid. De stof wordt in elk segment afgebroken, evenredig met de snelheid die er in elk segment is ( f i = K i ⋅ g i ). De uitwisseling tussen de segmenten is evenredig met de aanwezige hoeveelheid: f ij = k ij x j .
a. Maak een toestandsmodel. b. Wat zijn de toestanden en wat zijn de ingangen? 3
3. Via een infuus wordt vloeistof Wi [m / s ] met daarin opgelost mineraal de bloedbaan met volume
mi [kg / s] ingebracht in
VB [m 3 ] . Door transpiratie verdwijnt Wt [m 3 / s] vloeistof uit het bload. 3
De nieren verwijderen Wn [ m / s ] vloeistof inclusief het mineraal (bij transpiratie blijft het mineraal 3
in het bloed). De blaas was al gevuld met VBL [m ] vloeistof, voordat het infuus werd aangebracht. a. Stel een toestandmodel op. b. Bereken de concentratie van het mineraal, zowel in de bloedbaan als in de blaas. 30
c. Wat zijn de toestanden en wat zijn de ingangen? 4. Stel je bevat steeds 5 liter bloed. Je drinkt op
t = 0 50 gram alcohol, die direct in je bloed komt. De
−1
afbraaksnelheid is k = 2 [h ] (dus x& = − kx ). a. Bepaal het toestandsmodel voor de hoeveelheid alcohol in je bloed. b. Wanneer is er minder dan 5 gram alcohol over? Tip: als
y& = −cy en y(0) = B → y(t ) = Be − ct
5. Tijdens Carnaval wordt wel eens bier gedronken. Veronderstel dat de alcohol direct in het bloed komt. Het bloedvolume is ongeveer constant ( V liter). De alcohol wordt afgebroken met reactieconstante c [1 / s ] evenredig met massa alcohol ( − cm ). Op t = 0 is er geen alcohol in het bloed. 1 Glas bevat 18 gram alcohol, c = 0.01, T = 100 , V = 3 l . Gebruik Beschrijf de concentratie alcohol in je bloed als je:
δ (t ) voor drinken
a. 4 glazen tegelijk drinkt b. 4 glazen na elke T seconden drinkt. . Schets het verloop van de concentratie in beide gevallen.
4.8 Chemical models 1. Gegeven is de volgende reactievegelijking:
A+ B
k1 k2
C
t = 0 is er a, b resp. c [mol ] van elke stof. Noem x [mol ] het aandeel van A dat is omgezet of aangemaakt.
Op
x& . b. Bereken x SS als geldt dat k1 = k 2 = 0, a = 3, b = 2, c = 1 . c. Wat is de minimale en wat is de maximale waarde van x ? Het antwoord is onafhankelijk van k1 en k 2 ).
a. Bereken
k A⎯ ⎯→ B + C met a, b, c het aantal mol van A, B, C op t = 0 . Noem het aantal mol van C dat wordt gemaakt x . a. Bereken x& . b. Bereken x SS .
2. Gegeven de volgende reactievergelijking:
4.9 Linearization 1. Gegeven is het volgende niet-lineaire systeem:
x& (t ) = 2 x 2 (t ) − x 3 (t ) a. Bereken de werkpunten van dit systeem. b. Schets x& als functie van x . Verklaar aan de hand van deze grafiek of het lineaire model in een werkpunt stabiel is of niet. 2. Gegeven is het niet-lineaire systeem:
y& (t ) + y 2 (t ) + y (t ) − u (t ) ⋅ y (t ) = 0 a. Bereken een/de werkpunt(en) van dit systeem wanneer b. Bereken in elk werkpunt het lineaire model. c. Beredeneer of dat lineaire model stabiel is. 3. Gegeven is het volgende niet-lineaire model:
x& (t ) + x 2 (t ) = 4u 2 (t ) y& (t ) + y 3 (t ) = 4 x 2 (t ) 31
u (t ) = 2 + δu (t ) .
a. Bepaal de werkpunten van het niet-lineaire model als
u0 = 1 .
b. Bereken in elk werkpunt het gelineariseerde toestandsmodel. Twee vaten zijn onderling verbonden door een buis. Deze buis heeft een lineaire hydraulische weerstand: p = k12φ . In elk vat stroomt een flow
φi [m 3 / s ] . Vat 1 heeft ook een afvoer via een
buis. Deze afvoerbuis heeft een niet-lineaire hydraulische
p = k1φ 2 . De hydraulische capaciteit is van 5 beide vaten C [m / N ] . Kies C = 1, k12 = 1, k1 = 1, φ1 = φ 2 = 2 .
weerstand k1 :
a. Bereken het toestandsmodel. b. Bereken het werkpunt. c. Bereken een lineair model. 4. Gegeven is het volgende niet-lineaire systeem:
z&(t ) + z 2 (t ) + z (t ) = 2 a. Bereken een/de werkpunt(en) van dit systeem wanneer het systeem in rust is: b. Bereken in elk werkpunt het lineaire model. c. Beredeneer of dat lineaire model stabiel is.
z&(t ) = 0 .
5. Gegeven is het volgende niet-lineaire systeem:
y& (t ) + y 2 (t ) = 9 a. Bereken een/de werkpunt(en) van dit systeem wanneer het systeem in rust is: b. Bereken in elk werkpunt het lineaire model. c. Beredeneer of dat lineaire model stabiel is.
y& (t ) = 0 .
4.10 Transfer functions 1. Bereken, zo mogelijk met MATLAB, de analytische vorm van
Y (s) =
a.
y (t ) . Bijvoorbeeld: als
3 −t , dan y (t ) = 3e , t ≥ 0 . Er zijn geen complexe getallen in y (t ) toegestaan. s +1
Y (s) =
5 s 2 + 3s + 2 ( s + 1)( s + 2)( s + 3)
3s 2 + 4s + 7 s 2 + 3s + 2 s+3 c. Y ( s ) = 2 s + 2s + 1 2s + 6 d. Y ( s ) = 2 s + 4s + 8 b.
Y (s) =
2s 5 + 3s 4 + 4s 3 + 13 e. Y ( s ) = 6 s + 21s 5 + 175s 4 + 735s 3 + 1624s 2 + 1764s + 720 2. Bereken, zo mogelijk met MATLAB, het toestandsmodel van:
32
s+3 s + 3s + 2 5( s + 2) b. H ( s ) = 2 s + 4s + 3 s 2 + 2s + 7 c. H ( s ) = 4 s + 3s 3 + 4s 2 + 2s + 1 a.
H ( s) =
2
3. Bereken de differentiaalvergelijking tussen uitgang overdrachtsfuncties bij vraag 2. 4. Bereken
σ , ω, ωn , ζ
y (t ) en ingang u (t ) van de drie
van de volgende polen:
s1 = −3 b. s1, 2 = −1 ± j a.
c.
s1, 2 = 3 ± 5 j
5. Gegeven is een toestandsmodel. Bereken de overdrachtsfunctie.
⎧ ⎛ − 4 − 8⎞ ⎛1⎞ ⎟⎟ ⋅ x + ⎜⎜ ⎟⎟ ⋅ u ⎪ x& = ⎜⎜ a. ⎨ 0 ⎠ ⎝ 1 ⎝0⎠ ⎪ y = (2 3) ⋅ x ⎩ ⎧ ⎛ − 4 − 8⎞ ⎛0⎞ ⎟⎟ ⋅ x + ⎜⎜ ⎟⎟ ⋅ u ⎪ x& = ⎜⎜ b. ⎨ 1 ⎠ ⎝ 0 ⎝1⎠ ⎪ y = (2 3) ⋅ x ⎩
c.
⎧ ⎛− 4 0 ⎞ ⎛ 1⎞ ⎟⎟ ⋅ x + ⎜⎜ ⎟⎟ ⋅ u ⎪ x& = ⎜⎜ ⎨ ⎝ 0 − 8⎠ ⎝ 1⎠ ⎪ y = (2 3) ⋅ x ⎩
Wat zijn de polen van H (s ) en wat zijn de eigenwaarden van A bij a) , b) en c) ? Teken het simulatieschema in simulink van alle 3 systemen.
4.11 Frequency responses 1) Indien op een eerste-orde fysiologisch systeem ( y& + ay = bu ) een sinusvormig ingangssignaal
u (t ) = sin(ω1t ) wordt aangebracht, blijkt de uitgang een amplitude te hebben van 0.5 als ω1 = 10 [rad / s ] . Als de frequentie verdubbelt ( ω 2 = 2ω1 ), dan wordt de amplitude van y (t ) 0.25 . Bereken a en b .
33
2. In de figuur zijn u (t ) en y (t ) getekend. De 10
10 en die van y (t ) is 6 . T = 1 s en ∆T = 0.125 s . y (t ) is de uitgang van k met ingang het eerste orde proces H ( s ) = s+ p u (t ) . Bereken k en p .
amplitude van u (t ) is
u(t)
6 y(t)
t
∆Τ
T
4.12 Regulation 1. Een statisch schema voor de bloedsomloop is als volgt gegeven:
Bereken P, Q als
Rv = 10, Ra = 1, PH = 10 .
Bereken de invloed op ( P, Q ) als
P
Q Ra
Rv , Ra , PH
Q
elk met 10% toenemen (drie situaties dus, telkens verandert maar één parameter).
1/Rv
P
2. Voor de regeling van een hart- en longmachine wordt gebruik gemaakt van het volgende regelschema voor de bloeddruk: Deze druk P wordt gemeten en doorgegeven aan een regelaar C . Deze stuurt de bloedstroom φ zodat de druk P wordt beïnvloed. Er is een verstoring d aangebracht als gevolg van de invloed van de patiënt. Het proces H wordt beschreven met p& + ap = bφ . De regelaar C regelt volgens
φ = 50e .
d = 0 . Bereken P (s ) als functie van Pref (s) . b) Veronderstel d ≠ 0 , dus we hebben een D (s ) . Bereken P (s ) als functie van Pref (s) en D (s ) .
a) Veronderstel
c) Veronderstel Pref (t ) = Pr U s (t ) en
d (t ) = d ssU s (t ) . Bereken P ( s ) = P (∞ ) .
d) Stel a = 10, b = 10 . Bereken de invloed van e) Herhaal vraag d als
d ss en Pr op Pss .
a en b +10% veranderen (dus 1 situatie).
34
