PROYECTO FINAL: DISEÑO DE UN
CONTROL PARA EL PÉNDULO INVERTIDO
GIRATORIO QUANSER SRV02-ETS Y
ROTPEN.
Betancourt, Duvan. Casallas, Hamilton
Universidad de La Salle, Bogotá, Colombia
{dbetancourt49, hcasallas17}@[Link]
1. Introducción
Se realizara el diseño de un control por el método de espacio de estados para controlar la posición
angular del péndulo invertido de la marca Quanser. Para ello se hallan las características dinámicas
del sistema del péndulo, específicamente las ecuaciones que describen el movimiento del brazo
rotativo y el péndulo con respecto al voltaje en el servomotor SRV02.
2. Metodología.
2.1. Modelo matemático
Figura 1. Convenciones del péndulo invertido
El modelo de péndulo invertido giratorio se muestra en la Figura 1. El pivote del brazo giratorio está
unido al sistema SRV02 y se activa. El brazo tiene una longitud de Lr, un momento de inercia de Jr,
y su ángulo, θ, aumenta positivamente cuando gira en sentido anti horario (CCW). El servo (y, por
lo tanto, el brazo) debe girar en la dirección CCW cuando la tensión de control es positiva, es decir,
Vm 0 .
El enlace del péndulo está conectado al extremo del brazo giratorio. Tiene una longitud total de Lp
y su centro de masa es L p / 2 . El momento de inercia sobre su centro de masa es Jp. El ángulo del
péndulo invertido, α, es cero cuando está perfectamente en posición vertical en la posición vertical
y aumenta positivamente cuando se gira CCW.
2.2. Modelado de péndulo rotacional invertido
Se realiza el modelo lagrangiano, este se define así:
L K U (1)
En términos generales:
L(q(t ), q(t )) K (q(t ), q (t )) U (q (t )) (2)
Finalmente las ecuaciones de movimiento de Euler-Lagrange, para un sistema de n grados de
libertad está definida de la siguiente manera:
d L(q, q) L(q, q)
i (3)
dq qi qi
Para determinar el modelo, se siguen seis pasos:
a. Determinación ecuación energía cinética.
La energía cinética está definida por una componente traslacional y otra rotacional para el péndulo
y una componente rotacional para el brazo:
1 1
K mv 2 I 2 (4)
2 2
Donde m es la masa del cuerpo, v es la velocidad lineal, I el momento de inercia, es la velocidad
angular y K es la energía cinética. Así la energía cinética para el brazo es:
1
K arm J arm 2 (5)
2
La energía cinética para el péndulo es:
1 1
Kp J p 2 m p v p 2 (6)
2 2
La energía cinética total del sistema es:
1 1 1
KT J p 2 m p v p 2 J arm 2 (7)
2 2 2
b. Determinación de ecuación de energía potencial:
Este sistema sólo almacena energía potencial gravitacional en el péndulo:
U m p gh (8)
U m p gl p cos 1 (9)
c. Centro de masa del péndulo:
Debido a que el péndulo es un cuerpo rígido, la posición requerida es la de su centro de masa:
X x Lr cos l p sin sin (10)
X y Lr sin l p sin cos (11)
X z l p cos (12)
La posición del centro de masa es:
X CM X X XY XZ (13)
d. Velocidad lineal del CM
Se debe determinar la velocidad del centro de masa, para obtener la ecuación del lagrangiano; para
ello se deriva la posición respecto al tiempo:
d
X CM VCM (14)
dt
T
X CM VCM X X XY X Z (15)
Se deriva las ecuaciones (10), (11), (12):
X x Lr sin l p sin cos cos cos (16)
X y Lr cos l p cos cos sin sin (17)
X z l p sin (18)
Con la velocidad lineal del centro de masa se obtiene el término 𝑉𝐶𝑀2 mediante la expresión:
XX
VCM X X
2
XY X Z X Y X x 2 X y 2 X z 2 (19)
XZ
Se hallan por separado el cuadrado de las componentes de la velocidad del centro de masa:
X x 2 2 Lr 2 sin 2 2 2 Lr sin l p sin cos cos sin l p 2 sin cos cos sin 2 (20)
X y 2 2 Lr 2 cos 2 2 2 Lr cos l p cos cos sin sin l p 2 cos cos sin sin 2 (21)
X z 2 2l p 2 sin 2 (22)
Se reemplaza (20), (21) y (22) en (19) y se simplifica:
VCM 2 Lr 2 2 l p 2 2 sin 2 2 Lr l p m p cos (23)
e. Energía total se sistema:
La energía cinética total del sistema es:
KT K P Karm (24)
Para determinar la ecuación se debe reemplazar (6) y (7) en (24):
1 1 1
KT J arm 2 J p 2 m p vCM 2 (25)
2 2 2
Se sustituye (23) y (25) para obtener la energía cinética total, la energía potencial se determinó en
la ecuación (10):
J arm 2 J p 2 m p Lr 2 2 l p 2 2 sin 2 2 Lr l p m p cos
1 1 1
KT
2 2 2
KT J arm 2 J p 2 m p Lr 2 2 Lp 2 2 sin 2 Lr Lp m p cos
1 1 1 1
(26)
2 2 2 2
f. Determinación de las ecuaciones de Euler lagrange:
L KT UT (27)
De esta manera el lagrangiano se determina como:
J arm 2 J p 2 m p Lr 2 2 l p 2 2 sin 2 2 Lr l p m p cos m p gl p 1 cos (28)
1 1 1
L
2 2 2
2.3. Ecuaciones no lineales de movimiento (Quanser Inc., 2012).
La dinámica, se obtendrán usando la ecuación de Euler - Lagrange:
2 L L
Qi (29)
tqi qi
Las variables qi se llaman coordenadas generalizadas. Para este sistema son:
q t (t ) (t )
T
(30)
Donde, como se muestra en la Figura 1, θ (t) es el ángulo del brazo giratorio y α (t) es el ángulo del
péndulo invertido. Las velocidades correspondientes son:
t t
q t
T
(31)
t t
Con las coordenadas generalizadas definidas, las ecuaciones de Euler-Lagrange para el sistema de
péndulo rotatorio son:
2 L L
Q1 (32)
t
2 L L
Q1 (33)
t
Las fuerzas generalizadas Qi se utilizan para describir las fuerzas no conservativas (por ejemplo, la
fricción) aplicadas a un sistema con respecto a las coordenadas generalizadas. En este caso, la fuerza
generalizada que actúa sobre el brazo giratorio y el péndulo es:
Q1 Br (34)
Q2 Bp (35)
Las ecuaciones de Euler-Lagrange son un método sistemático para encontrar las ecuaciones de
movimiento, es decir, las EOM de un sistema. Después de pasar por este proceso, las ecuaciones de
movimiento no lineales para el SRV02 péndulo invertido rotativo son:
1 1 1 1 2
m p Lr m p Lp m p Lp cos J r m p Lp sin cos m p Lp Lr sin Br
2 2 2 2 2
(36)
4 4 2 2
1 1 1 1
m p Lp 2 sin cos J p m p Lp 2 m p Lp 2 cos cos m p Lp 2 g sin B p (37)
2 4 4 2
El par motor se aplica a la base del brazo giratorio (es decir, en el engranaje de carga) como se
describe por la ecuación:
g K gm kt Vm K g km
(38)
Rm
Si se observa ambas ecuaciones coincides con la forma típica de una MOE para un solo cuerpo:
Jx bx gx 1
Donde x es una posición angular, J es el momento de inercia, b es el amortiguamiento, g es la función
gravitacional y τ1 es el par aplicado es un valor escalar.
2.4. Linealización
Modelo de espacio de estado lineal
Las ecuaciones lineales de espacio de estado son:
x Ax Bu (39)
y Cx Du
Donde x es el estado, u es la entrada de control, A, B, C y D son matrices de espacio de estado. Para
el sistema de péndulo rotatorio, se definen el estado y la salida:
xT (40)
y x1 x2
Se procede con la linealización de la primera ecuación del péndulo giratorio invertido no
lineal (36), definiendo las condiciones iniciales de las variables en cero, así:
z T 0 0 0 0 0 0 (41)
De esta manera la ecuación (36) se linealizada para y da (las otras derivadas dan 0):
f ( z ) 1 1
m p L2r m p L2p m p L2p cos(0) J r m p L2r J r (42)
z z 0 4 4
f ( z ) 1 1
m p Lp Lr cos(0) J r m p Lp Lr (43)
z z 0 2 2
Al aplicar la ecuación de linealización completa:
f ( z ) f ( z ) f ( z ) f ( z ) f ( z ) f ( z )
flin( z )) f ( z0 ) (43)
z z 0 z z 0 z z 0 z z 0 z z 0 z z 0
Se obtiene:
m L p
2
r J r m p Lp Lr B p
1
2
(44)
Se procede con la linealización de la ecuación del péndulo giratorio invertido no lineal (37),
definiendo las condiciones iniciales de las variables en cero, así:
z T 0 0 0 0 0 0 (45)
De esta manera la ecuación (36) se linealizada para , y da (las otras derivadas dan cero):
f ( z ) 1
m p Lp Lr (46)
z z 0 2
f ( z ) 1
J p m p L2p (47)
z z 0 4
f ( z ) 1
m p Lp g (48)
z z 0 2
Al aplicar la ecuación de linealización completa, se obtiene:
1 1 1
m p Lp Lr J p m p L2p m p L p g B p (49)
2 4 2
Se ajustan las dos ecuaciones lineales de movimiento encontradas en la forma de matriz
que se muestra en la Ecuación (50), forma de MOE para dos cuerpos:
d11 d12 c c g
d q 11 12 q 1 1 (50)
21 d22 c21 c22 g 2 2
1
m p Lr J r m p Lr Lp 0
2
0
2 Br
(51)
Bp m p Lp g 0
1
1 m L L 0
J p m p L2p
1
2 p r p 2
4
Se resuelven los términos de aceleración en las ecuaciones de movimiento, se hace usando
la inversa de la matriz:
1
a b
1 1 d b
A c a (52)
c d det( A)
El determinante de A:
1
J T m p L2r J r J p m p L2p m 2p L2p L2r
1
(53)
4 4
Resolviendo:
1 1
J p m p L2p m p Lp Lr Br
1 4 2
1
m L g B
(54)
JT 1 m L L 2
J r m p Lr 2 p p
p
2 p p r
Expandir la ecuación y recopilar términos semejantes:
1 1 2 1 1 2 2 1 2
J p 4 m p Lp Br 2 m p Lp Lr B p 4 m p Lp Lr g J p 4 m p L p (55)
JT
1 1
m p Lp Lr Br J r m p L2r B p m p Lp g J r m p L2r m p L p Lr (56)
1 1
JT 2 2 2
2.5. Espacio de estados
Finalmente se encuentra el espacio de estados lineal del sistema de péndulo invertido rotatorio.
Donde:
x1
x2
x3
x4
De esta manera la construcción de las matrices A, B, C Y D del espacio de estados queda determinada
así:
0 0 1 0
0 0 0 1
1 2 2 1 1
A 0 m p Lp g J P mP L2p Br mP L p Lr (57)
4 4 2
1
J r mP L2r B p
1 1
0 mP Lp g J r mP L2r mP L p Lr Br
2 4 2
0
0
1 1 1 0 0 0 0
B J P mP Lp C
2
D
JT
4 0 1 0 0 0
1m L L
2 P p r
2.6. Diseño del control
Figura 2. Esquema de control por realimentación de estados.
Una vez obtenido el espacio de estados, se procede con el diseño del controlador, para el Quanser
se hace uso del método de retroalimentación de estados. Para realizar la simulación en primer lugar
se definen los parámetros según las especificaciones del fabricante, como se muestra en la siguiente
tabla:
Figura 3. Especificaciones técnicas del Quanser (Quanser Inc., 2012).
Posteriormente se definen las matrices del espacio de estados, tal como se observa en la ecuación
57, donde además se agregan parámetros referentes al funcionamiento del servomotor SRV02
presente. Se debe tener en cuenta que tanto la masa del péndulo mp, el momento de inercia del
péndulo respecto a sus centro de masa Jp,cm y el momento de inercia del brazo respecto a su centro
de masa Jarm,cm, cambian dependiendo si se agregan o no masas extra al péndulo. El espacio de
estados para el sistema sin agregar masas es el siguiente:
0 0 1 0
0 0 0 1
A (57)
0 81.3358 -46.0337 -0.9312
0 121.9664 -46.0337 -1.3964
0
1 0 1 0 0 0 0
B C D
J T 83.6100 0 1 0 0 0
80.4395
Una vez definido el espacio de estados se procede con la definición de los parámetros de diseño. En
cuanto al polinomio de diseño al que se igualara la función de transferencia del diagrama completo,
se determinan sus parámetros de la siguiente manera:
R1 600, R2 200. Mp 0.1, Tss 0.35
2
Mp
Ln
100 , Wn
4
(58)
Mp
2
Tss
Ln
2
100
Quedando determina la función de transferencia diseño con la ecuación 1.3:
Sal s 157.6
(59)
Pd s s 4 822.9s 3 1.384e05s 2 2.869e06s 1.892e07
Luego se discretiza tanto el espacio de estados del sistema, como la función de diseño con un tiempo
de discretización de 0.002 segundos; ahora se ha determinado la matriz A y B discretas y los polos
de diseño discreto, por ende se puede usar la función del programa Matlab llamada place, de
manera que se pueden determinar las ganancias K, que serán usadas tanto en la simulación, como
en la implementación del control.
Se debe recordar que el vector de estados está determinado así:
x (60)
El proceso de diseño se realiza de nuevo pero esta vez con varias masas, que se van agregando en
siete instantes de tiempo, dichas piezas sin impresas en ABS: el acrilonitrilo butadieno estireno a 40
% de densidad. Se debe tener en cuenta que cuando se agregan nuevas masas estas cambian el
momento de inercia del péndulo y la masa del sistema péndulo-masas insecto. Ahora se procede
con la determinación del momento de inercia con masas adicionales:
Primero se determina el momento de inercia del péndulo que en si es una barrilla uniforme:
Figura 4. Barra rígida.
L/ 2
M
L/ 2 M r3 M L3 L3
I r dr
2
L/ 2 L L 3 L/ 2
3 L 8 8
1
I CM ML2
12
El momento de inercia de una barra con una masa (Ver Figura 4), en el caso del péndulo se
determina así:
1
I P ,CM M p L p 2 M i ha 2 ( 61 )
12
Donde Mi es la masa agregada y ha es la altura desde el pivote al centro de la masa. En las siguientes
imágenes se muestran los tipos de masas agregadas al péndulo.
Figura 5. Masa cilíndrica con agujero interno (3 gramos).
Figura 6. Argolla de cobre (11,6 gramos).
Control Switcheado:
Las siete ganancias halladas se muestran a continuación:
K 7 [-1.7756 21.2865 -1.9848 2.4605] : Instante uno, sin masa adicional
K 6 [-2.3965 27.8874 -2.4758 3.2895] : Instante dos, masa 3 gramos
K5 = [-2.9881 34.1733 -2.9429 4.0912]: Instante tres, masa 6 gramos
K 4 =[-3.3469 38.0304 -3.2264 4.6045] : Instante cuatro, masa 9 gramos
K3 [-3.6748 41.5369 -3.4849 5.0675] : Instante cinco, masa 12 gramos
K 2 [-5.4204 57.3966 -4.7874 6.7214]: Instante seis, masa 15 gramos
K1 [-7.3643 75.2383 -6.2398 8.6417]: Instante siete, masa 26.6 gramos
Ahora se decide implementar un switcheado, para ello se hace uso de las ganancias anteriores; por
medio de ciclos if – else se switchea cada una de las ganancias para diferentes rangos de par Angulo
del péndulo ( Angulo) los rangos comienzan en 7 grados y aumentan de dos en dos, exceptuando
el rango siete que posee un rango desde 17 a 20. La implementación del código de control
switcheado se muestra como pseudocódigo en la figura 9.
2.7. Pseudocódigos
Modeling
Limpiar consola , Borrar todo
Variables Mp(cambiantes), Lp ,lp , Jp(cambiantes), Dp , marm , r , Lr , Jarm _ cm , Dr ,Kenc, g,Jt
Inicio
Definición Matriz A R n, B n 1, C 1 n , D 0
Se agrupan las matricesA, B ,C , D a espacio de estados
Se discretiza el espacio de estados, con tustin a Ts= 0.01
Def overshoot Mp 1%
2
Mp
Ln
Def de coeficiente de amortiguamiento C 1 100
2
Mp
2 Ln
100
4
Def de frecuencia natural Wn
Ts ee
Wn 2
Def de tf de diseño H 2
s 2CWn Wn2 s 10 s 11
Se discretiza la tf de diseño : Hd
Se obtienen los los polos de diseño de la tf Hd : polos
Deter min ar la ganancia K place(Ad ,Cd , polos)
Variables x , [ 1.5;0; 1;0 ] y, u( señal entrada ), k(index)
Hacer
u(k) ref Kx(k)
y(k 1 ) Cd x(k) Dd u(k) rand
x(k 1 ) Ad x(k) Bd u(k)
Hasta que k tamaño de t 1
Grafica posicion pendulo Vs tiempo(t)
Grafica velocidad pendulo Vs tiempo(t)
Grafica posicion motor Vs tiempo(t)
Grafica velocidad motor Vs tiempo(t)
Grafica voltaje Vs tiempo(t)
Fin del programa
Figura 7. Pseudocódigo para simulación del control diseñado.
Para la simulación del control se hace uso de la herramienta Matlab, teniendo en cuenta, cada uno
de los vectores K obtenidos, en siete instantes, donde cada instante vario así: primero se evalúa y
determina un control para el péndulo si masa adicional, luego con una masa adicional de 3 gramos
a 25 cm del pivote, de ahí en adelante se van agregando masas de 3 gramos cada una, la distancia
ha se disminuye en un centímetro, ya que este es el ancho de la arandela plástica, finalmente se
agrega una masa de cobre de 11,6 gramos a 17 cm respecto al pivote.
En cuanto al código de implementación, pseudocódigo de este se presenta a continuación:
Deprontosi
Limpiar consola , Borrar todo
Variables Mp(Instante 1 ), Lp ,lp , Jp(Instante 1 ), Dp , marm , r , Lr , Jarm _ cm , Dr ,Kenc,g, Jt
Inicio
Definición Matriz A R n, B n 1, C 1 n , D 0
Se agrupan las matricesA, B ,C , D a espacio de estados
Se discretiza el espacio de estados, con tustin a Ts= 0.01
Def overshoot Mp 1%
2
Mp
Ln
Def de coeficiente de amortiguamiento C 1 100
2
Mp
2 Ln
100
4
Def de frecuencia natural Wn
Ts ee
Wn 2
Def de tf de diseño H 2
s 2CWn Wn2 s 10 s 11
Se discretiza la tf de diseño : Hd
Se obtienen los los polos de diseño de la tf Hd : polos
Deter min ar la ganancia K place(Ad ,Cd , polos)
try
Definición de tipo de placa , reloj , frecuencia , muestreo por segundo ,...
canal de lectura , canal de escritura , amplitud , ganancia , periodo ,t
Apertura de comunicación
Establecer encoder : canal 0,1
Identificar placa
Establecer rango de salida
Crear , niciar y set mascara de datos para lectura encoderdel motor
Crear , iniciar y set mascara de datos para lectura encoderdel pendulo
Definición var iable pmotor , pmotor _ ant , vmotor , vmotor _ ant , vpendulo ,...
vpendulo _ ant , ppendulo _ arriba , ppendulo _ ant , ppendulo , almacenamiento
Definición del vector de estados x
Mientras
Actualizar estados : anterior actual
Lectura enco der del motor
Conversión de conteos a radianes
Calcular mod ulo a 2 radianes : ppednulo _ arriba
Se deriva la posición : velocidad del motor
Lectura enco der del pendulo
Conversión de conteos a radianes
Se deriva la posición : velocidad del pendulo
Se det er min a el vector de espacio de estados
si posición del pendulo 0.2618 rad
voltaje Kx : acción de control
Se colocan saturadores de voltaje inverso
Escritura de voltaje en canal ana log ico
sin o
voltaje es cero
Escritura de voltaje en canal ana log ico
Fin
aumenta i , guardar el tiempo , posición , velocidad de motor y pendulo
Hasta que Esc no se oprima
Cerrar y borrar mascara de lecturade motor y pendulo
Apagar motor y cerrar comunicación
Grafica posicion pendulo Vs tiempo(t)
Grafica velocidad pendulo Vs tiempo(t)
Grafica posicion motor Vs tiempo(t)
Grafica velocidad motor Vs tiempo(t)
Grafica voltaje Vs tiempo(t)
catch error
Fin del programa
Figura 8. Pseudocódigo para implementación del control diseñado sin agregar masa.
Deprontosi
Limpiar consola, Borrar todo
Se colocan los vectores de ganancia K1 , K 2 , K 3 , K 4 , K 5 , K 6 , K 7
Inicio
try
Definición de tipo de placa, reloj , frecuencia , muestreo por segundo,...
canal de lectura, canal de escritura , amplitud , ganancia, periodo, t
Apertura de comunicación
Establecer encoder : canal 0,1
Identificar placa
Establecer rango de salida
Crear , niciar y set mascara de datos para lectura encoderdel motor
Crear , iniciar y set mascara de datos para lectura encoderdel pendulo
Definición var iable pmotor , pmotor _ ant , vmotor , vmotor _ ant , vpendulo ,...
vpendulo _ ant , ppendulo _ arriba , ppendulo _ ant , ppendulo, almacenamiento
Definición del vector de estados x
Mientras
Actualizar estados : anterior actual
Lectura enco der del motor
Conversión de conteos a radianes
Calcular mod ulo a 2 radianes : ppednulo _ arriba
Se deriva la posición : velocidad del motor
Lectura enco der del pendulo
Conversión de conteos a radianes
Se deriva la posición : velocidad del pendulo
Se det er min a el vector de espacio de estados
si posición del pendulo 0.349066 rad y posición del pendulo 0.296706 rad
voltaje K 7 x : acción de control
sin o posición del pendulo 0.296706 rad y posición del pendulo 0.261799 rad
voltaje K 6 x : acción de control
sin o posición del pendulo 0.261799 rad y posición del pendulo 0.226893 rad
voltaje K 5 x : acción de control
sin o posición del pendulo 0.226893 rad y posición del pendulo 0.191986 rad
voltaje K 4 x : acción de control
sin o posición del pendulo 0.191986 rad y posición del pendulo 0.15708 rad
voltaje K 3 x : acción de control
sin o posición del pendulo 0.15708 rad y posición del pendulo 0.122173 rad
voltaje K 2 x : acción de control
sin o posición del pendulo 0.122173 rad
voltaje K1 x : acción de control
Se colocan saturadores de voltaje inverso
Escritura de voltaje en canal ana log ico
sin o
voltaje es cero
Escritura de voltaje en canal ana log ico
Fin
aumenta i, guardar el tiempo, posición, velocidad de motor y pendulo
Hasta que Esc no se oprima
Cerrar y borrar mascara de lecturade motor y pendulo
Apagar motor y cerrar comunicación
Grafica posicion pendulo Vs tiempo (t)
Grafica velocidad pendulo Vs tiempo (t)
Grafica posicion motor Vs tiempo (t)
Grafica velocidad motor Vs tiempo(t)
Grafica voltaje Vs tiempo (t)
catch error
Fin del programa
3.1. Simulación
Se realiza básicamente simulaciones para el péndulo sin masa, con masas de 3 y 11,6 gramos:
Figura 9. Simulación de péndulo invertido sin masa adicional (Instante 1).
Figura 10. Simulación de péndulo invertido con masa de 3 gramos (Instante 2).
Figura 11. Simulación de péndulo invertido con masa de 15 gramos (Instante 6).
Figura 12. Simulación de péndulo invertido si con masa de 26,6 gramos (Instante 7).
3.2. Implementación
Para la implementación se realiza el siguiente montaje:
Figura 13. Montaje del Quanser.
Figura 14. Graficas de velocidad, posición y voltaje en la implementación del control.
Figura 15. Prueba de funcionamiento de control.
4. Conclusiones
Mientras se realizaban pruebas de implementación se determina que uno de los cuatro
estados, exactamente es el estado de posición del motor es inestable.
Para obtener la velocidad tanto del péndulo como del motor, se hace uso de un filtro
discreto, dicho filtro debe ser discretizado con el mismo tiempo que se usa para discretizar
el espacio de estados y la función de transferencia de diseño, el cual es 2 milésimas de
segundo.
El resultado del control switcheado implementado es más robusto que el control sencillo,
esto se debe al cambio entre constantes de realimentación, a medida que se cambia la
posición angular del péndulo.
5. Bibliografía
Quanser Inc. (2012). Inverted Pendulum Experiment. USER MANUAL. Markham, Ontario, Canada:
Quanser Inc.
Quanser Inc. (2012). Inverted Pendulum Experiment for MatLab®/simulink® users. Instructor
Workbook. Markham, Ontario, Canada: Quanser Inc.