UNIVERSIDAD NACIONAL DE

SAN AGUSTÍN
FACULTAD DE INGENIERÍA DE
PRODUCCIÓN Y SERVICIOS

Escuela Profesional de

U Ingeniería Electrónica

Curso:

N Teoría de control automático 1

Tema:
Guías de laboratorio

S Elaborado por:
Ing William Vladimir Mullisaca Atamari
Ing Oscar Salazar Alarcón

A Ing Juan Carlos Cuadros Machuca

Dr Daniel Yanyachi Aco Cardenas

Arequipa – Perú

2019
 TEMARIO DE LAS GUÍAS DE LABORATORIO

Laboratorio 1: Introducción Al Matlab Para Sistemas de control lineal.

Solucionario de laboratorio 1

Laboratorio 2: Modelado De Sistemas Físicos Para Sistemas De Control Lineal.

Solucionario de laboratorio 2

Laboratorio 3: Simulación De Sistemas Numéricos De Sistemas De Control Lineal.

Solucionario de laboratorio 3

Laboratorio 4: La Respuesta Temporal De Los Sistemas De Control Lineal.

Solucionario de laboratorio 4

Laboratorio 5: Respuesta En Frecuencia De Los Sistemas De Control Lineal.

Solucionario de laboratorio 5
 LABORATORIO 1

I OBJETIVO

I.1. Emplear MATLAB en la solución de problemas matemáticos que involucren a la Transformada

de Laplace y la Transformada Inversa de Laplace.

II PROCEDIMIENTO
II.1. Repetir y ejercitar los siguientes comandos en Matlab. Definición de una constante:

a=1
b= [1 2]

Escribiendo números complejos:

a=2+i
b=-5-3*i

Expresión booleana:

a==1

Vector constante:

v=[1 2 3 4 5]
ó v=1:5

Matriz constante:

A= [2 2 3
007
5 9 -1]
ó
A= [2 2 3;0 0 7;5 9 -1]

Podemos formar matrices usando operaciones con objetos definidos anteriormente:

a=1; b=2;

Observe que, si colocamos punto y coma al final de la expresión, el resultado no es mostrado en la

pantalla, lo que puede ser conveniente en algunas situaciones.

A= [a+b pi 3
b^2 0 atan(a)
5 sin(b) -1]

Podemos formar matrices y vectores de zeros:

B=zeros()
Matriz de zeros con 2 las y 3 columnas:

B=zeros(2,3)

Matriz de zeros con las dimensiones de la matriz A:

A= [2 2 3;0 0 7;5 9 -1]; B=zeros(A)

De modo semejante, podemos formar matrices y vectores de unos: Matriz de unos con 2 filas y 3
columnas:

C=ones(2,3)

Matrices diagonales:

Matriz diagonal con los elementos de la diagonal principal yendo de 1 a 5:

D=diag(1:5)

Extrayendo los elementos de la diagonal principal:

A= [1 2 3
456
7 8 9]
B=diag(A)

Formando una matriz diagonal con los elementos de la diagonal principal de una matriz:

C=diag(diag(A))

Operaciones: Matriz identidad:

A=diag(ones(1,3))
ó
A=eye(3)

Suma de matrices (recuerde las matrices deben tener la misma dimensión):

B=A+A

Sumar 1 a todos los elementos de una matriz:

C=B+1

Multiplicación de matrices:

A= [1 2 3;4 5 6;7 8 9]
C= [1 2 0;0 0 1;0 2 3]
D=A*C

Multiplicación elemento a elemento:

A= [1 0 0;0 2 3;5 0 4]
B= [2 0 0;0 2 2;0 0 3]
C=A.*B

Extracción de la fila 2:

a=C (2, :)

Extracción de la columna 3:

b=C (:,3)

Traza de una matriz:

A= [1 2 3; 4 5 6; 7 8 9]
t=trace(A)

Rank (rango) de una matriz:

r=rank(A)

Matriz transpuesta:

B=A'

Inversa de una matriz:

A= [0 1; -2 -3]
B=inv(A)
A*B

Determinante de una matriz:

d=det(A)

Polinomios:

Polinomio p1 con raíces en 0 y -1:

v= [0 -1]
p1=poly(v)

Polinomio p2 con coeficientes 1 y 2 y 1:

p2=poly([1 2 1])

Cálculo de raíces:

p=roots(p1)

Autovalores y autovectores (eigenvalues and eigenvectors):

 Autovalores:

A= [0 1; -2 -3]
r=eig(A) % r vector de autovalores ó
[V,D]=eig(A) % produce un matriz diagonal D de autovalores
% y una matriz completa V cuyas columnas son
% sus correspondientes autovectores. Así (A*V=V*D)

Funciones:

function[y]=mifuncion(x) % definimos un archivo mifuncion.m

if x<0 % donde estará nuestra 'y=f(x)'

y=x^2

else

y=sin(x*(pi/180))

end

Una vez definida la función, podemos evaluarla en el punto x=30 (grados)

y=mifuncion(30)

Otra función (mifuncion2): la función y = x + x2 + sin((2 ∗ pi) ∗ x) en Matlab.

Plotear la función entre -2 y 3:

Cree un vector variando de -2 a 3 con paso de 0.5

x=-2:0.5:3;
Calculando la función
y=mifuncion2(x)
Ploteando el resultado:
plot(x,y)

Creando nuevas ventanas gráficas (figure(1), figure(2) ... figure(n)):

figure(1) % abre nueva ventana de gráfico

plot(x,y,'r') % la curva es de color rojo
figure(2) % abre nueva ventana de gráfico
plot(x,y,'b:') % la curva es de color azul y entrepunteada
III Ejercicios:

Implemente en matlab la siguiente función, luego plotee:

Implemente en matlab la siguiente función, luego plotee:

Haga un .m file que ayude a encontrar el mínimo de f (x) = x3 − 2x − 5, dentro del intervalo
(0,2)

Construya una señal escalón unitario, de 0 a 50 segundos, con step inicial en 25 s. El paso
deberá ser de 0.5s. Plotee el resultado

Construya una señal delta de dirac, plotee el resultado.

IV MATERIAL Y EQUIPO

IV.1. Una PC con SO Windows XP y MATLAB

V OBSERVACIONES Y CONCLUSIONES

V.1. Haga sus observaciones y emita al menos cinco conclusiones del trabajo realizado
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
 SOLUCIÓN LABORATORIO 1

Ejercicios:
1.- Implemente en matlab la siguiente función, luego plotee:

1
𝑦 = 𝑓(𝑥) =
𝑥2 −1

Ejecutando en Matlab

Paso 1: Escribimos en Editor:

clear,clc
x=-10:0.1:10; %Definimos los valores de x
y=1./(x.^2-1); %Definimos la función

plot(x,y) %Grafica de f(x)

grid %Ponemos cuadricula y nombres
a ejes para mejor visualización
title('y=1/(x^2-1)')
xlabel('Valores x')
ylabel('valores y')

Fig. 1.1 Ventana Editor

Paso 2: Ejecutamos Run en Editor (ver Fig. 1.1)

 Fig.1.2 Ploteo después de ejecutar Run

2.- Implemente en Matlab la siguiente función, luego plotee:

𝑥1 + 𝑥2 𝑖𝑓 𝑥1 > 0, 𝑥2 > 0
𝑦 = 𝑓(𝑥1 , 𝑥2 ) = {
√𝑥12 + 𝑥22 𝑒𝑛 𝑙𝑜𝑠 𝑑𝑒𝑚á𝑠 𝑐𝑎𝑠𝑜𝑠

Ejecutando en Matlab

Paso1: Escribimos en Editor

clear,clc
x=-10:0.2:10; %Definimos x
y=-10:0.2:10; %Definimos y
[X,Y]=meshgrid(x,y); %Usaremos surf por ende
necesitamos matrices cuadradas X,Y necesitamos matrices cuadradas X,Y
%meshgrid(x,y) repetirá los elementos de x en filas, y en columnas guardadas en [X,Y] respectivamente

Z=zeros(1,101); %preasignación de la cantidad

máxima de espacio requerido para la
formación.

for i=1:length(x) %Iniciamos un bucle para operar

todas las combinaciones posibles

for j=1:length(y) %x1,x2 son de la misma

magnitud 101x101
a=X(i,j); %Guardamos cada valor
ubicado de X en a
b=Y(i,j); %Guardamos cada valor
ubicado de Y en b

if a>0 && b>0 %Hacemos la primera

comparación si a>0 y b>0

Z(i,j)=a+b; %Colocamos la función

 condicional
else
Z(i,j)=(a^2+b^2)^(1/2); %Colocamos la
función para
los demas casos
end %Terminamos los bucles

end

end

surf(X,Y,Z) %Ploteamos

title('Z = F(X,Y)') %Colocamos título y nombres

a ejes para mejor visualización

xlabel('Valores X'),ylabel('Valores Y'),zlabel('Valores Z')

Fig. 2.1 Programa para poner graficar función por parte f(x1,x2)

Paso 2: ejecutar Run

Fig. 2.2 Ploteo de f(x1,x2) por partes

3.- Haga un .m file que ayude a encontrar el mínimo de f (x) = x3 − 2x − 5, dentro del intervalo
(0,2)

Ejecutando en Matlab

Paso1: Escribimos en Editor

clear,clc

for x=0:0.0001:2 %Definimos x de 0 a 2 con

intervalo muy pequeño para mayor exactitud

y=x.^3-x.*2-5; %Definimos la función y

if x==0 %Guardamos el primer valor de

x,y como mínimos
xm=x;
ym=y;
end

if ym>y %Comparamos el valor mínimo

con el valor y
ym=y;
xm=x;
end
end

str=['El mínimos valor se encuentra en x=',num2str(xm),', y=',num2str(ym),''];

disp(str) %Mostramos la respuesta en la

ventana de comandos

Fig. 3.1 Programa para encontrar el mínimo valor de un polinomio

Paso 2: Ejecutar Run

 Fig. 3.2 Obtención del valor minino en ventana comando

4.- Construya una señal escalón unitario, de 0 a 50 segundos, con step inicial en 25 s. El paso
deberá ser de 0.5s. Plotee el resultado

Ejecutando en Matlab

Paso1: Escribir el programa en Editor:

clear,clc

t=0:0.5:50; %Definimos dominio t

T=25; %Definimos la posición del escalón

y=stepfun(t,T); %Usamos stepfun como función

para escalón

plot(t,y),grid; %Graficamos
title('Escalón Unitario')
xlabel('Tiempo (t)')
ylabel('Funcion f(t)')

axis([0 50 -0.5 1.5]) %Usamos axis para apreciar

mejor el gráfico

Fig. 4.1 Programa para graficar función escalón unitario en t=25s

Paso 2: ejecutar Run

Fig. 4.2 Ploteo del escalón unitario con salto en t=25s

5.- Construya una señal peine de dirac, plotee el resultado.

Ejecutando en Matlab:
Paso 1: escribimos en Editor

clear,clc;

x=-10:1:10 % valores de x, con salto de uno

y=zeros(1,length(x)); %matrix fila cero

s=0; %valor de inicio de s

for i=1:length(x);
pi=dirac(x-x(i)); %pi=valor de Inf en la
posición i

if s==0;
s=pi;
end

s=s+pi; %suma de todas las matrices

formadas

end

% cdv=Cambio de valor Inf por valor real para poder graficar:

cdv=s==Inf; % encontrando Inf

y(cdv)=1; % valor de y (puede variar de
0 a infinito)
stem(x,y,'^','fill') % plotea
title('Peine de Dirac')
xlabel('Tiempo (t)')
ylabel('Funcion f(t)')

Fig. 5.1 programa para el Peine de Dirac

Paso 2: Ejecutar Run

Fig. 5.2 Tren de pulso o Peine de Dirac con intervalo de un 1s

6.- OBSERVACIONES Y CONCLUSIONES

 Para la creación de programas (Funcion1, ...funcion5).se ha usado la ventana Editor de Matlab

 Las funciones f(x) de una variable; para la obtención de su grafica se usa el comando plot(x,y) y
para funciones f(x,y) de dos variables el comando surf(X,Y,Z) previo uso del comando
 [X,Y]=meshgrid(x,y) para obtener datos matriciales que darán los valores de Z mediante comando
for-end y if-end para funciones por partes como es y=f(x1´x2)
 El comando for-end permite hacer varias iteraciones, lo cual permite encontrar valores máximos o
mínimos de polinomios algebraicos
 Stepfun(T,T0) es un comando para obtener la función escalón unitario con salto en T0 cualquiera.
 El uso adecuado de for-end y if-end nos permite a partir de una función como es Dirac(comando de
Matlab que da la función delta dirac)obtener una nueva función (peine de dirac)
 Para la gráfica de pulsos (como Dirac, previo cambio de valor de Inf con valor real) en un punto
dado se usa el comando stem(x,y).
 LABORATORIO 2

I OBJETIVO

I.1. Aprovechar el software MATLAB en la solución de ecuaciones diferenciales, en la solución de

problemas matemáticos en el área de control.

II ASPECTOS TEÓRICOS

 Solución de ecuaciones diferenciales, modelamiento de sistemas físicos.

III MATERIAL Y EQUIPO

III.1. Una PC con SO Windows XP y MATLAB

IV PROCEDIMIENTO

1. Dada la siguiente ecuación diferencial:

Hallar la solución por el método que usted considere apropiado, también debe graficar la
solución.

2. Ecuaciones diferenciales (Repetir la siguiente actividad)

Primero, implementamos la función en edo1.m

Function [ydot]=edol(y)
ydot=(1-y)/2;

Segundo, cree el archivo euler_ode1.m

t(1)=0; % Instante inicial

tf=10; % Instante final
 y(1)=0; % Condición inicial
ye(1)=0; % Valor inicial de la solución exacta
%
% Paso de integración (experimente alterar el paso):
h=0.5;
% Cálculo de número de pasos):
n=round(tf/h);
% Integración numérica usando el método de Euler:
% Comando for:
for i=1:n
% Vector de tiempo:
t(i+1)=t(i)+h;
% Vector ydot (derivada en el tiempo de y):
ydot(i)=edol(y(i));
% solución numérica:
y(i+1)=y(i)+h*ydot(i);
% Solución exacta:
ye(i+1)=1-exp(-t(i+1)/2);
% Finaliza comando for:
end
% Determinación del último término del vector ydot:
ydot(n+1)=edol(y(n+1));
% Ploteando la solución exacta 'ye' y la solución numérica 'y', ambos versus el vector
plot(t,y,t,ye,'r:');
% Colocando una leyenda en la parte superior derecha de la figura:
legend('solución exacta','solución numérica - Euler')
% Colocando título en la parte superior y etiquetas en las coordenadas title('Comparación
entre la solución exacta y la solución numérica')
xlabel('Tiempo (s)')
ylabel('Amplitud ()')

Para ejecutar la simulación, en el workspace digite euler_ode1

Obs: Los archivos .m deben ser guardados en el directorio de trabajo, que generalmente es:
C:\MATLAB6p5\work
Comentario: Uno de los métodos numéricos más usados es el método Runge Kutta de cuarta
orden (RK4), que brinda soluciones más exactas en comparación con el método Euler.

3. Dado el siguiente circuito eléctrico, hallar la ecuación diferencial que lo representa:

Figura 1: Circuito RC

4. De la figura 1 anterior hallar su representación en diagrama de bloques y su función de

transferencia.

5. De la figura 1 hallar su representación mediante diagramas de flujo de señal y obtenga

 su función de transferencia mediante la regla de mason.

6. De la figura 1, obtenga la solución de la ecuación diferencial usando la herramienta

MATLAB.

7. Ejercicio: reducir el diagrama de bloques de la figura 2 y expresar la función de

transferencia entre C y R en función de (G1, G2 y H)

Figura 2: Diagrama de Bloques.

V OBSERVACIONES Y CONCLUSIONES

V.1. Haga sus observaciones y emita al menos cinco conclusiones del trabajo realizado
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
 SOLUCIÓN LABORATORIO 2

1. Dada la siguiente ecuación diferencial:

[𝟏 − 𝒚(𝒕)]
𝒚̇ =
𝟐
Hallar la solución por el método que usted considere apropiado, también debe graficar la
solución.

Aplicando transformada de Laplace:

[1/S − Y]
SY =
2
2SY = 1/S − Y
1/S 1 1
Y= = −
(2S + 1) 𝑆 𝑆 + 1/2

Transformada de Laplace inversa:

𝐲(𝐭) = 𝟏 − ℮(−𝐭/𝟐)

Fig. 1 función y (t)

2. Ecuaciones diferenciales (Repetir la siguiente actividad)

[𝟏−𝒚(𝒕)]
• Determine la solución ecuación diferencial 𝒚̇ = 𝟐 usando el método de Euler.
• 𝑷𝒍𝒐𝒕𝒆𝒂𝒓 el grafico de 𝒚(𝒕) en función de t, y el grafico de la solución exacta 𝒚𝒆 (𝒕)en
función de t.

Instante inicial: 𝒕 = 𝟎
Instante final: 𝒕𝒇 = 𝟏𝟎
Solución exacta de la ecuación diferencial: 𝒚𝒆 (𝒕) = 𝟏 − ℮(−𝒕/𝟐)

Resolviendo por medio de 𝒎𝒂𝒕𝒍𝒂𝒃

a.- Creando en la ventana editor:
 edo1.m
function[𝑦𝑑𝑜𝑡] = 𝑒𝑑𝑜1(𝑦)
𝑦𝑑𝑜𝑡 = (1 − 𝑦)/2

b.- Creando en la ventana editor: euler_edo1.m

c.- Ejecutando la simulación; presionando: 𝐑𝐮𝐧 𝐞𝐮𝐥𝐞𝐫_𝐞𝐝𝐨𝟏. 𝐦

comparación entre la solución exacta y la solución numérica

1

0.9
solución numérica - euler
0.8
solución exacta
0.7

0.6
Amplitud ()

0.5

0.4

0.3

0.2

0.1

0
0 1 2 3 4 5 6 7 8 9 10
Tiempo (s)

3. Dado el siguiente circuito eléctrico, hallar la ecuación diferencial que lo representa:

 Figura 1: Circuito RC
Malla i

R(i − 𝑖2 ) = e

Malla i2

1 𝑡
R(𝑖2 − i) + ∫ 𝑖 dt = 0
𝐶 0 2

4. De la figura 1 anterior hallar su representación en diagrama de bloques y su función de

transferencia.

Usando transformada de Laplace en las ecuaciones anteriores:

Malla I
R(I − 𝐼2 ) = E
E
I = 𝐼2 +
R

Malla I2
𝐼2
R(𝐼2 − I) + =0
𝐶𝑆
1
(R + ) 𝐼 = RI
𝐶𝑆 2
1
(1 + )𝐼 = I
𝑅𝐶𝑆 2
1
𝐼2 = I
1
(1 + )
𝑅𝐶𝑆
 𝑅𝐶𝑆
𝐼2 = I
𝑅𝐶𝑆 + 1

Resultado:

RCS
𝐼2 RCS +1 ∗1
=
𝐸 RCS R
1−
RCS + 1
𝐼2
= 𝐶𝑆
𝐸

5. De la figura 1 hallar su representación mediante diagramas de flujo de señal y obtenga

su función de transferencia mediante la regla de Mason.

𝑉2 𝑃𝑘 ∆𝑘
=∑
𝑉1 ∆
1 RCS
𝑃1 = ∗
R RCS + 1
RCS
𝐿1 =
RCS + 1
RCS
∆= 1 −
RCS + 1
∆1 = 1
1 RCS
𝐼2 R ∗ RCS + 1
=
𝐸 1 − RCS
RCS + 1
𝐼2
= 𝐶𝑆
𝐸
6. De la figura 1, obtenga la solución de la ecuación diferencial usando la herramienta
MATLAB.

7. Ejercicio: reducir el diagrama de bloques de la figura 2 y expresar la función de

transferencia entre C y R en función de (G1, G2 y H)

Figura 2: Diagrama de Bloques.

Resolución por medio del algebra de bloques:

Moviendo el punto de ramificación 1 a la derecha del bloque G2(S) OBTENEMOS:
RESULTADO:
C(S) G2(S) + G1(S)
=
R(S) 1 + G2(S)H(S)

CONCLUSIONES
 El uso del método de Euler es un método práctico para la comprensión de otros métodos matemático
en el campo de ecuaciones diferenciales.
 El uso adecuado del método de diagramas de bloque algebraico permite tener un resultado correcto
para sistemas de control.
 LABORATORIO 3

I OBJETIVO
I.1. Emplear MATLAB en la solución de problemas matemáticos en el área de control.

II ASPECTOS TEÓRICOS

 Dados en clase de teoría.

III MATERIAL Y EQUIPO

III.1. Una PC con SO Windows XP y MATLAB

IV PROCEDIMIENTO

1. Dada la ecuación diferencial, resuelta en el laboratorio 2.

Y recordando como hallamos su solución numérica en Matlab.

“…Primero, implementamos la función en edo1.m

Function [ydot]=edo1(y)
ydot=(1-y) /2;

Segundo, cree el archivo euler_ode1.m

t(1)=0; % Instante inicial

tf=10; % Instante final
y(1)=0; % Condición inicial
ye(1)=0; % Valor inicial de la solución exacta
%
% Paso de integración (experimente alterar el paso):
h=0.5;
% Cálculo de número de pasos):
n=round(tf/h);
% Integración numérica usando el método de Euler:
% Comando for:
for i=1:n
% Vector de tiempo:
t(i+1)=t(i)+h;
% Vector ydot (derivada en el tiempo de y):
ydot(i)=edol(y(i));
% solución numérica:
y(i+1)=y(i)+h*ydot(i);
 % Solución exacta:
ye(i+1)=1-exp(-t(i+1)/2);
% Finaliza comando for:
end
% Determinación del ultimo termino del vector ydot:
ydot(n+1)=edol(y(n+1));
% Ploteando la solución exacta 'ye' y la solución numérica 'y', ambos versus el vector
plot(t,y,t,ye,'r:');
% Colocando una leyenda en la parte superior derecha de la figura:
legend('solución exacta','solución numérica - Euler')
% Colocando título en la parte superior y etiquetas en las coordenadas
title('Comparación entre la solución exacta y la solución numérica')
xlabel('Tiempo (s)')
ylabel('Amplitud ()')…”

Para ejecutar la simulación, en el workspace digite euler_ode1

2. Del ítem 1, repita las instrucciones, le ayudara a entender como introducir funciones de
transferencia (LTI) en matlab:

Definimos nuestro sistema expresado por la F.T. H(s) = 1/(s + 1)

h=tf(1,[1 1])

Definimos nuestro sistema G expresado en espacio de estados:

x˙ = Ax + Bu
y = Cx + Du
g=ss([0 1;-5 -2],[0;3],[0 1],0)

Obtención del numerador y denominador de la F.T. a partir de las matrices A, B, C y D (matrices

del modelo en espacio de estados). [num,den]=ss2tf([0 1;-5 -2],[0;3],[0 1],0)

Principales gráficas para facilitar el análisis de sistemas LTI.

g = tf([16],[1 9 16]) % definimos nuestro sistema subplot(2,2,1)

bode(g) % grafico de bode subplot(2,2,2)
pzmap(g) % grafico de polos y zeros subplot(2,2,3)
step(g) % respuesta escalón unitario subplot(2,2,4)
nyquist(g) % grafico de Nyquist

3. Definir que es un espacio de estados, escriba un ejemplo, y cual es su aplicación en

teoría de control automático.

4. Ejercicio: realice la simulación numérica de la ecuación diferencial descrita en el item

1. Está vez implemente el algoritmo RK4, en vez del método Euler. Haga una comparación
entre las soluciones: Exacta, Euler y RK4.

5. Ejercicio: repetir 1, usando la ecuación diferencial (no-lineal) siguiente:

 y˙ = 2.5 [y(1 − y)] y(0) = 0.9 Condición inicial h = 0.85 Paso tf = 600s Instante final

6. Ejercicio: el propósito de estos ejercicios es familiarizarse con las funciones de

transferencia, su representación y respuesta. Ud. podrá ver la respuesta temporal. Por
tanto, para las siguientes funciones de transferencia grafique (i) el diagrama de polos y
zeros, (ii) la respuesta escalón unitario, (iii) la respuesta a la señal rampa:

V OBSERVACIONES Y CONCLUSIONES

V.1. Haga sus observaciones y emita al menos cinco conclusiones del trabajo realizado
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
 SOLUCIÓN LABORATORIO 3

1. Dada la ecuación diferencial, resuelta en el laboratorio 2.

Y recordando como hallamos su solución numérica en Matlab.

“…Primero, implementamos la función en edo1.m

function [ydot]=edo1(y)
ydot=(1-y)/2;

Segundo, cree el archivo euler_ode1.m

t(1)=0; % Instante inicial

tf=10; % Instante final
y(1)=0; % Condición inicial
ye(1)=0; % Valor inicial de la solución exacta
% Paso de integración (experimente alterar el paso):
h=0.5;
% Cálculo de número de pasos):
n=round(tf/h);
% Integración numérica usando el método de Euler:
% Comando for:
for i=1:n
% Vector de tiempo:
t(i+1)=t(i)+h;
% Vector ydot (derivada en el tiempo de y):
ydot(i)=edo1(y(i));
% solución numérica:
y(i+1)=y(i)+h*ydot(i);
% Solución exacta:
ye(i+1)=1-exp(-t(i+1)/2);
% Finaliza comando for:
end
% Determinación del ultimo termino del vector ydot:
ydot(n+1)=edo1(y(n+1));
% Ploteando la solución exacta 'ye' y la solución numérica 'y', ambos versus el vector
plot(t,y,t,ye,'r:');
% Colocando una leyenda en la parte superior derecha de la figura:
legend('solución exacta','solución numérica - Euler')
% Colocando título en la parte superior y etiquetas en las coordenadas title('Comparación entre
la solución exacta y la solución numérica')
xlabel('Tiempo (s)')
ylabel('Amplitud ()')

Para ejecutar la simulación, en el workspace digite euler_ode1

2. Del ítem 1, repita las instrucciones, le ayudara a entender como introducir funciones de
transferencia (LTI) en matlab:

Definimos nuestro sistema expresado por la F.T. H(s) = 1/(s + 1)

h=tf(1,[1 1])

Definimos nuestro sistema G expresado en espacio de estados:

x˙ = Ax + Bu
y = Cx + Du
g=ss([0 1;-5 -2],[0;3],[0 1],0)

Obtención del numerador y denominador de la F.T. a partir de las matrices A, B, C y D (matrices

del modelo en espacio de estados).
[num,den]=ss2tf([0 1;-5 -2],[0;3],[0 1],0)

Principales gráficas para facilitar el análisis de sistemas LTI.

g = tf([16],[1 9 16]) % definimos nuestro sistema subplot(2,2,1)

bode(g) % grafico de bode subplot(2,2,2)
pzmap(g) % grafico de polos y zeros subplot(2,2,3)
step(g) % respuesta escalón unitario subplot(2,2,4)
nyquist(g) % grafico de Nyquist

3. Definir que es un espacio de estados, escriba un ejemplo, y cuál es su aplicación en

teoría de control automático.

 Espacio de Estados: Un modelo matemático de un sistema físico descrito mediante un entradas,

salidas y variables de estado relacionadas por ecuaciones diferenciales que se combinan en
 una ecuación diferencial matricial de primer orden.

La representación en espacio de estados puede hallarse través de la función de tranferencia.

Las variables son expresadas como vectores (entradas, salidas y estados) y las ecuaciones
algebraicas se escriben en forma matricial (esto último solo cuando el sistema dinámico es lineal
e invariante en el tiempo).

 Variables de estado: Las variables de estado son el subconjunto más pequeño de variables de
un sistema que pueden representar su estado dinámico completo en un determinado instante. El
número mínimo de variables de estado necesarias para representar un sistema dado, n, es
normalmente igual al orden de la ecuación diferencial que define al sistema.

 Ejemplo de Espacio de Estados con función de transferencia

 Ejemplo de Espacio de Estados con ayuda de Matlab.

  Los espacios de estados, en lo que es Control, se utiliza tanto para observar el
comportamiento de un sistema físico dinámico transformado a ecuaciones diferenciales como
para controlar la respuesta de acuerdo al valor de las entradas.

4. Ejercicio: realice la simulación numérica de la ecuación diferencial descrita en el item 1.

Está vez implemente el algoritmo RK4, en vez del método Euler. Haga una comparación
entre las soluciones: Exacta, Euler y RK4.

 Para resolver consultamos en internet sobre el método RK4

 Agregaremos el siguiente código a euler_ode1.m para reciclar.

 Para separar los métodos usamos nuestro valor inicial con una nueva variable:
yRK4(1)=0; % Valor inicial para metodo RK4

 Ahora agregamos el método al bucle ‘for’, donde se desarrollan los demás métodos, teniendo en
cuenta que nuestra ecuación edo1.m no se encuentra en función de x, por lo tanto eliminamos
esta variable de nuestro método:
 % Solucion RK4:
k1=edo1(yRK4(i));
k2=edo1(yRK4(i)+k1/2);
k3=edo1(yRK4(i)+k2/2);
k4=edo1(yRK4(i)+k3);
yRK4(i+1)=yRK4(i)+h*(k1+2*k2+2*k3+k4)/6;

 Ahora para poder analizar los 3 métodos juntos modificamos la siguiente línea agregando la
visualización del método RK4
% Ploteando la solución exacta 'ye', la solución numérica 'y' y la solución Runge-Kutta 4°orden
'yRK4' en versus
plot(t,y,t,ye,'r:',t,yRK4,'b:');
% Colocando una leyenda en la parte superior derecha de la figura:
legend('solución exacta','solución numérica - Euler','solución RK4')

 Ahora guardamos nuestro nuevo código bajo el nombre euler_rk4_ode1.m y ejecutamos.

  El método RK4 es menos exacta que la solución de Euler y más lejana a la exacta; debido a que
esta funciona con una ED de 4to orden y con 2 variables, teniendo en cuenta que nuestra edo1.m
es una ED de 1er orden y 1 variable.

5. Ejercicio: repetir 1, usando la ecuación diferencial (no-lineal) siguiente:

y˙ = 2.5 [y(1 − y)]

y(0) = 0.9 Condición inicial
h = 0.85 Paso
tf = 600s Instante final

 Creamos un nuevo script con la ecuación diferencial no lineal ed_nl.m

function [ydot]=ed_nl(y)
ydot=2.5*(y*(1-y));

 Modificamos los parámetros iniciales de euler_ode1.m bajo el nombre euler_nl.m

t(1)=0; % Instante inicial

tf=600; % Instante final
y(1)=0.9; % Condición inicial
ye(1)=0.9; % Valor inicial de la solución exacta
% Paso de integración (experimente alterar el paso):
h=0.85;
% Cálculo de número de pasos:
n=round(tf/h);
% Integración numérica usando el método de Euler:
% Comando for:
for i=1:n
% Vector de tiempo:
t(i+1)=t(i)+h;
% Vector ydot (derivada en el tiempo de y):
ydot(i)=ed_nl(y(i));
% solución numérica:
 y(i+1)=y(i)+h*ydot(i);
% Solución exacta:
ye(i+1)=1-exp(-t(i+1)/2);
% Finaliza comando for:
end
% Determinación del ultimo termino del vector ydot:
ydot(n+1)=ed_nl(y(n+1));
% Ploteando la solución exacta 'ye' y la solución numérica 'y', ambos versus el vector
plot(t,y,t,ye,'r:'),axis([0 600 0.2 1.2]);
% Colocando una leyenda en la parte superior derecha de la figura:
legend('solución exacta','solución numérica - Euler')
% Colocando título en la parte superior y etiquetas en las coordenadas title('Comparación entre la
solución exacta y la solución numérica')
xlabel('Tiempo (s)')
ylabel('Amplitud ()')

 Ejecutamos el programa

 Los métodos no son exactos ya que están diseñados para ecuaciones de primer grado
ordinarios y lineales.

6. Ejercicio: el propósito de estos ejercicios es familiarizarse con las funciones de

transferencia, su representación y respuesta. Ud. podrá ver la respuesta temporal. Por
tanto, para las siguientes funciones de transferencia grafique (i) el diagrama de polos y
zeros, (ii) la respuesta escalón unitario, (iii) la respuesta a la señal rampa:
a) Sistema de primer orden

%a%
num=1;
den=[1 1];
g=tf(num,den);

t=0:0.1:10;
ramp=t;
r=lsim(g,ramp,t);
ts=ones(1,length(t));

subplot(2,1,1);plot(t,step(g,t),t,ts,'b:',t,r,'r',t,t,'r:'),axis([0 10 0 4])
grid on
legend('Respuesta al escalon','Escalon','Respuesta a la rampa','Rampa')
xlabel('Tiempo(s)'),ylabel('Amplitud')
subplot(2,1,2),pzmap(g),
b) Sistema de segundo orden con 2 polos

%b%
num=2;
den1=[1 1];
den2=[1 2];
den=conv(den1,den2);
g=tf(num,den);

t=0:0.1:5;
ramp=t;
r=lsim(g,ramp,t);
ts=ones(1,length(t));

subplot(2,1,1);plot(t,step(g,t),t,ts,'b:',t,r,'r',t,t,'r:'),axis([0 5 0 3]), grid on

legend('Respuesta al escalon','Escalon','Respuesta a la rampa','Rampa')
xlabel('Tiempo(s)'),ylabel('Amplitud')
subplot(2,1,2),pzmap(g),
c) Sistemas de primer orden con zero en semiplano derecho e izquierdo.
%c%
num1=[1 -2]; num2=[1 2];
den=[2 1];
g1=tf(num1,den); g2=tf(num2,den);

t=0:0.1:10;
ramp=t;
r1=lsim(g1,ramp,t);
ts=ones(1,length(t));
r2=lsim(g2,ramp,t);

subplot(2,2,1);plot(t,step(g1,t),t,ts,'b:',t,r1,'r',t,t,'r:'),axis([0 10 -10 2]),grid on

legend('Respuesta al escalon','Escalon','Respuesta a la
rampa','Rampa'),xlabel('Tiempo(s)'),ylabel('Amplitud'),title('H1')
subplot(2,2,3),pzmap(g1)

subplot(2,2,2);plot(t,step(g2,t),t,ts,'b:',t,r2,'r',t,t,'r:'),axis([0 10 0 10])
grid on
legend('Respuesta al escalon','Escalon','Respuesta a la
rampa','Rampa'),xlabel('Tiempo(s)'),ylabel('Amplitud'),title('H2')
subplot(2,2,4),pzmap(g2)
d) Sistema de segundo orden con Wn=5 y ζ=0.2

%d%
num=25;
den=[1 2 25];
g=tf(num,den);

t=0:0.01:5;
ramp=t;
r=lsim(g,ramp,t);
ts=ones(1,length(t));

subplot(2,1,1);plot(t,step(g,t),t,ts,'b:',t,r,'r',t,t,'r:'),axis([0 5 0 3])
grid on
legend('Respuesta al escalon','Escalon','Respuesta a la rampa','Rampa')
xlabel('Tiempo(s)'),ylabel('Amplitud')
subplot(2,1,2),pzmap(g)
e) Sistema de inciso d) con un polo en el origen

%e%
num=[25 0];
den=[1 2 25];
g=tf(num,den);

t=0:0.01:5;
ramp=t;
r=lsim(g,ramp,t);
ts=ones(1,length(t));

subplot(2,1,1);plot(t,step(g,t),t,ts,'b:',t,r,'r',t,t,'r:')
grid on
legend('Respuesta al escalon','Escalon','Respuesta a la rampa','Rampa')
xlabel('Tiempo(s)'),ylabel('Amplitud')
subplot(2,1,2),pzmap(g),
f) Sistema de inciso e) con cero en -5 rad/s

%f%
num1=5;
num2=[1 5];
num=conv(num1,num2);
den=[1 2 25];

g=tf(num,den);

t=0:0.01:5;
ramp=t;
r=lsim(g,ramp,t);
ts=ones(1,length(t));

subplot(2,1,1);plot(t,step(g,t),t,ts,'b:',t,r,'r',t,t,'r:'),axis([0 5 0 5])
grid on
legend('Respuesta al escalon','Escalon','Respuesta a la rampa','Rampa')
xlabel('Tiempo(s)'),ylabel('Amplitud')
subplot(2,1,2),pzmap(g),
g) Sistema de inciso f) con un tercer polo a baja y alta frecuencia
%g%
num1=25; num2=125;
denx=[1 2 25];
den1=conv([1 1],denx); den2=conv([1 5],denx);
g1=tf(num1,den1); g2=tf(num2,den2);

t=0:0.01:5;
ramp=t;
r1=lsim(g1,ramp,t);
ts=ones(1,length(t));
r2=lsim(g2,ramp,t);

subplot(2,2,1);plot(t,step(g1,t),t,ts,'b:',t,r1,'r',t,t,'r:'),grid on
legend('Respuesta al escalon','Escalon','Respuesta a la rampa','Rampa')
xlabel('Tiempo(s)'),ylabel('Amplitud'),title('H1')
subplot(2,2,3),pzmap(g1)

subplot(2,2,2);plot(t,step(g2,t),t,ts,'b:',t,r2,'r',t,t,'r:')
grid on
legend('Respuesta al escalon','Escalon','Respuesta a la rampa','Rampa')
xlabel('Tiempo(s)'),ylabel('Amplitud'),title('H2')
subplot(2,2,4),pzmap(g2)
VI OBSERVACIONES Y CONCLUSIONES
 Observaciones:

- La respuesta en rampa no estaba incluida en Matlab por lo que se ideó un código para poder
plotear.

 Conclusiones:

- Se ha podido observar el cambio de una FT respecto a lo que son sus polos y ceros, y también
se ha podido estimar visualmente cual es su efecto al agregar al sistema.
- Respecto a lo anterior dicho se puede decir que la presencia de ceros puede modificar el valor
final que se tomara al analizar en régimen permanente.
- La presencia de polos inestabiliza la forma de onda en régimen transitorio al aplicar frecuencias
grandes.
- Se ha comprendido mas lo que son espacios de estados.
- El método de Euler y método exacto descrita al inicio en el inciso 1. Se emplean para dar
solución ecuaciones diferenciales ordinarias de primer orden (EDO 1er orden) lineales, por lo
cual, al aplicar una ecuación no lineal, hubo fallos en la solución.
 LABORATORIO 4

I OBJETIVO
I.1. Emplear MATLAB en la solución de problemas matemáticos en el área de control.

II ASPECTOS TEÓRICOS

 Dados en clase de teoría.

III MATERIAL Y EQUIPO

III.1. Una PC con SO Windows XP y MATLAB

IV PROCEDIMIENTO

Controlador PID: Repita las siguientes instrucciones y tome atención en cada uno de los resultados
obtenidos. Le permitirá conocer la acción proporcional (P), Integral (I) y derivativa (D), de los
controladores. Estructura mayormente usada en la industria.

El modelo dinámico de un motor DC está expresado por la siguiente función de transferencia:

Los valores numéricos de los coeficientes del modelo están en la Tabla 1.

Las especificaciones de diseño del controlador para este sistema, considerando una entrada escalón
unitario de 1 rad=s, son:

 tiempo de acomodación menor que 0.04 s.

 overshoot menor que 16 %.
 ess = 0 error en régimen permanente (steady state error)
La función de transferencia de nuestro controlador PID es de la forma:

Para resolver este problema, crearemos un _le pidmotor.m donde introduciremos lo siguiente:

Teniendo ya la estructura del sistema controlado, procedemos a sintonizar (tunning) los parámetros del
PID (Ecuación (4))

a) Acción Proporcional (P)

Observe. ¿Cuál es el efecto de la acción proporcional?

b) Acción Proporcional Integral (PI)

Incremente el valor de KI, primero a 50, luego a 200.

Observe. ¿Cuál es el efecto de la acción integral?
En los casos anteriores, el tiempo de acomodación es muy largo. Para reducirlo y obtener una
respuesta más rápida, hacemos las modificaciones KP = 17 y KI = 200.

c) Acción Proporcional Integral Derivativo (PID)

En el último caso, obtuvimos una respuesta rápida, pero el overshoot incremento demasiado. Por
tanto, para reducir este efecto, introduciremos la acción derivativa.

Con los datos KP = 17, KI = 200 y KD = 0:15 obtenemos:

Observe. ¿Cuál es el efecto de la acción derivativa?

Finalmente hemos alcanzado las especificaciones de diseño de controlador. El resultado muestra
un overshoot por debajo de 16 %, un tiempo de acomodación inferior a 40 ms y no hay error en
régimen permanente.
V OBSERVACIONES Y CONCLUSIONES
V.1. Haga sus observaciones y emita al menos cinco conclusiones del trabajo realizado
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
 SOLUCIÓN LABORATORIO 4

Controlador PID: Repita las siguientes instrucciones y tome atención en cada uno de los resultados
obtenidos. Le permitirá conocer la acción proporcional (P), Integral (I) y derivativa (D), de los
controladores.

Estructura mayormente usada en la industria.

El modelo dinámico de un motor DC está expresado por la siguiente función de transferencia:

Figura 4.1: Modelo dinámico de un motor DC

𝒔(𝑱𝒔 + 𝒃)𝛉(𝒔) = 𝒌𝑰(𝒔) (𝟏)

(𝑳𝒔 + 𝑹)𝑰(𝒔) = 𝑽 − 𝒌𝒔𝜽(𝒔) (𝟐)
𝜽̇ 𝒌
𝑮(𝒔) = = (𝟑)
𝑽 (𝑱𝒔 + 𝒃)(𝑳𝒔 + 𝑹) + 𝒌𝟐

Los valores numéricos de los coeficientes del modelo están en la Tabla 1.

Las especificaciones de diseño del controlador para este sistema, considerando una entrada escalón
unitario de 1rad=s, son:

• tiempo de acomodación menor que 0.04 s.

• overshoot menor que 16 %.
• ess = 0 error en régimen permanente (steady state error)

Tabla 4.1: Coeficiente para el modelo dinámico del motor DC

constante unidades descripción
J 0.01 kgm2/s2 Momento de inercia del motor
b 𝟎. 𝟏 𝑵𝒎/𝒔 Coeficiente de amortecimiento del
𝒌 = 𝒌𝒆 = 𝒌𝒕 𝟎. 𝟎𝟏 𝑵𝒎/𝑨𝒎𝒑 motor
R 𝟏𝛀 Constante de la fuerza electromotriz
L 0.5 H Resistencia eléctrica del motor
V ? Inductancia eléctrica del motor
? Voltaje de entrada
Ángulo o posición de salida

La función de transferencia de nuestro controlador PID es de la forma:

𝑲𝑰 𝑲𝑫 𝒔𝟐 + 𝑲𝑷 𝒔 + 𝑲𝑰
𝑷𝑰𝑫(𝒔) = 𝑲𝒑 + + 𝑲𝑫 𝒔 = (𝟒)
𝒔 𝒔

Para resolver este problema, crearemos un file 𝐩𝐢𝐝𝐦𝐨𝐭𝐨𝐫. 𝐦 donde introduciremos lo siguiente:
Teniendo ya la estructura del sistema controlado, procedemos a sintonizar (tunning) los parámetros
del PID (Ecuación (4))

a) Acción Proporcional (P)

clear all
h=pidmotor(1.79,0,0)
t=0:0.001:0.3; %tiempo de la simulación
step(h,t)
xlabel('Tiempo (s)')
ylabel('Posición (rad)')
grid
%Info de la respuesta escalon%
info=stepinfo(h);
tr=info.RiseTime;tp=info.PeakTime;ts=info.SettlingTime;Mp=info.Overshoot;
str=['tr=',num2str(tr),' ts=',num2str(ts),' tp=',num2str(tp),' Mp=',num2str(Mp)];
title(str)

Fig. 1
Observe. ¿Cuál es el efecto de la acción proporcional?

 El efecto proporcional de Kp es:

o Producir una overshoot (Mp) más alto.
o Aumentar el factor de amortiguamiento.
o El tiempo de asentamiento (ts) no varía considerablemente.
o No varia el error en estado estacionario. ess=0.
 El resultado es un régimen transitorio más inestable u oscilante, los efectos son proporcionales
al valor de Kp respecto a 1.
 Adicionalmente podemos decir que las especificaciones de diseño establecidos ts<0.04s y
Mp<16% no se cumplen. El error ess=0, se cumple.

b) Acción Proporcional Integral (PI)

clear all
h=pidmotor(1.79,20,0);
t=0:0.001:0.3; %tiempo de la simulación
step(h,t)
xlabel('Tiempo (s)')
ylabel('Posición (rad)')
grid
%info de la respuesta%
info=stepinfo(h);
tr=info.RiseTime;tp=info.PeakTime;ts=info.SettlingTime;Mp=info.Overshoot;
str=['tr=',num2str(tr),' ts=',num2str(ts),' tp=',num2str(tp),' Mp=',num2str(Mp)];
title(str)

Fig. 2

Incremente el valor de KI, primero a 50, luego a 200.

clear all
h=pidmotor(1.79,50,0);
t=0:0.001:0.3; %tiempo de la simulación
step(h,t)
xlabel('Tiempo (s)')
ylabel('Posición (rad)')
grid
%info de la respuesta%
info=stepinfo(h);
tr=info.RiseTime;tp=info.PeakTime;ts=info.SettlingTime;Mp=info.Overshoot;
str=['tr=',num2str(tr),' ts=',num2str(ts),' tp=',num2str(tp),' Mp=',num2str(Mp)];
title(str)

Fig. 3
clear all
h=pidmotor(1.79,50,0);
t=0:0.001:0.3; %tiempo de la simulación
step(h,t)
xlabel('Tiempo (s)')
ylabel('Posición (rad)')
grid
%info de la respuesta%
info=stepinfo(h);
tr=info.RiseTime;tp=info.PeakTime;ts=info.SettlingTime;Mp=info.Overshoot;
str=['tr=',num2str(tr),' ts=',num2str(ts),' tp=',num2str(tp),' Mp=',num2str(Mp)];
title(str)

Fig. 4
Observe. ¿Cuál es el efecto de la acción integral?
 El efecto de la acción integral (Ki):
o Aumenta el overshoot (Mp) hasta un 100% como máximo; y luego de que ‘ts’ se vuelva
infinito el error en estado estacionario es infinito.
o Reduce el factor de amortiguamiento y puede volverse negativo, obteniendo así un
tiempo de asentamiento (ts) mucho más alto y desde cierto punto llegando hasta
infinito.
o Varía el error en estado estacionario cambiando de ess=0 a, 1 e infinito
o Todo esto proporcional al valor Ki.

En los casos anteriores, el tiempo de acomodación es muy largo. Para reducirlo y obtener una
respuesta más rápida, hacemos las modificaciones KP = 17 y KI = 200.

clear all
h=pidmotor(17,200,0);
t=0:0.001:0.3; %tiempo de la simulación
step(h,t)
xlabel('Tiempo (s)')
ylabel('Posición (rad)')
grid
%info de la respuesta%
info=stepinfo(h);
tr=info.RiseTime;tp=info.PeakTime;ts=info.SettlingTime;Mp=info.Overshoot;
str=['tr=',num2str(tr),' ts=',num2str(ts),' tp=',num2str(tp),' Mp=',num2str(Mp)];
title(str)

Fig. 5
 Comentario:
- Aun no se cumplen las especificaciones de diseño Mp<16% y ts<0.04s
- El error ess=0, se cumple.

c) Acción Proporcional Integral Derivativo (PID)

 En el último caso, obtuvimos una respuesta rápida, pero el overshoot incremento demasiado. Por
tanto, para reducir este efecto, introduciremos la acción derivativa. Con los datos KP = 17, KI = 200
y KD = 0.15 obtenemos:

clear all
h=pidmotor(17,200,0.15);
t=0:0.001:0.3; %tiempo de la simulación
step(h,t)
xlabel('Tiempo (s)')
ylabel('Posición (rad)')
grid
%info de la respuesta%
info=stepinfo(h);
tr=info.RiseTime;tp=info.PeakTime;ts=info.SettlingTime;Mp=info.Overshoot;
str=['tr=',num2str(tr),' ts=',num2str(ts),' tp=',num2str(tp),' Mp=',num2str(Mp)];
title(str)

Fig. 6

Observe. ¿Cuál es el efecto de la acción derivativa?

 El efecto de la acción derivativa (K) es:

o Reduce drásticamente el overshoot (Mp).
o Aumenta drásticamente el factor de amortiguamiento.
o Reduce el tiempo de asentamiento (ts).
o No varia el error en estado estacionario. ess=0

Finalmente hemos alcanzad las especificaciones de diseño de controlador. El resultado muestra un

overshoot (Mp) por debajo de 16%, un tiempo de acomodación (ts) inferior a 40 ms y no hay error en
régimen permanente.

OBSERVACIONES Y CONCLUSIONES

Haga sus observaciones y emita al menos cinco conclusiones del trabajo realizado

 Los tipos de control estudiados tienen sus propias características y son adecuadas para
determinada aplicación, la aplicación de un controlador va a depender básicamente del tipo de
señal de entrada al sistema de control, en este caso un escalón unitario.
 Con el uso del MATLAB, podemos encontrar la respuesta temporal de un sistema de control
basándonos en su función de transferencia en LA o LC.
 Al variar los valores de 𝜍 y 𝑤𝑛 con el controlador PID, podemos hacer que el sistema de control
responda más rápidamente, o con menor rebase, menor oscilación, lo cual nos permite adecuar
el controlador a los requerimientos dado.
 Los controles PID reúnen las características de los tres tipos de control precedentes, por tanto, su
respuesta temporal es más idónea.
 El control proporcional, nos ayuda a definir el comportamiento en régimen transitorio.
 El control integral, nos permite definir el comportamiento tanto transitorio como en permanente
pudiendo hasta cierto punto inestabilizarlo, integrando el sistema.
 El control derivativo hace exactamente lo opuesto a la integral de manera mucho más rápida.
Según lo experimentado, estabiliza la respuesta transitoria y permanente, derivando el sistema.
 LABORATORIO 5

I OBJETIVO
I.1. Emplear MATLAB en la solución de problemas matemáticos en el área de control.

II ASPECTOS TEÓRICOS

 Dados en clase de teoría.

III MATERIAL Y EQUIPO

III.1. Una PC con SO Windows XP y MATLAB

IV PROCEDIMIENTO

1. Para fines didácticos, iremos analizar la respuesta en frecuencia de un sistema de suspensión de

automóvil.

%//////////////////////////////////////////////
% Sistema de Suspensión //
% /////////////////////////////////////////////
clear % borra las variables antiguas
% /////////////////////////////////////////////
% Parámetros
m=250; % masa suspensa (kg)
k=10000; % rigidez del resorte (N/m)
b=316.23; % amortecimiento (Ns/m)
% Matrices del sistema
A=[0 1;-k/m -b/m];
B=[-1;b/m];
C=[-k/m -b/m];
D=[b/m];
 % Sistema en el espacio de estados
G=ss(A,B,C,D);
% Respuesta en frecuencia
% Diagrama de Bode del sistema (observe las escalas de los graficos):
bode(G)
% La escala horizontal es logarítmica, y muestra la frecuencia en rad/s.
% La escala vertical del grafico superior es lineal, y muestra la ganancia en db.
% La escala vertical del grafico inferior es lineal y muestra la fase en grados.
% Matlab asigna valores por defecto a nuestro grafico de Bode, el Matlab
% mostro que la amplitud de la ganancia esta entre -0db e +40db, la
% fase está entre -90 e 90 grados, y que la frecuencia esta entre 10 y
% 1000 rad/s. Vamos construir un vector w donde el logaritmo de la
% distancia entre un valor y otro es constante, comenzando en 10
% elevado a (-3), terminando en 10 elevado a (4), con 2000 puntos, y
% asumiendo que la unidad es rad/s:
w=logspace(-3,4,2000);
% Determinamos los vetores de magnitud y fase,
[mag,phase]=bode(G,w);
% mag y fase son arrays de 3 dimensiones, esto quiere decir que la
% función bode trabaja con sistemas multivariable (MIMO)
% Entonces, para nuestro caso, debemos solo capturar el vector correspondiente
mag=mag(1,:);
phase=phase(1,:);
% Ploteando el Bode de ganancia
subplot(2,1,1)
semilogx(w,20*log10(mag)) % semilogx realiza un plot semilogaritmico
axis([10^-3 10^4 -70 60]) % definimos los márgenes de la ventana
grid % rejilla
% Ploteando el Bode de fase
subplot(2,1,2)
semilogx(w,phase) % semilogx realiza un plot semilogaritmico
axis([10^-3 10^4 -100 100]) % definimos los márgenes de la ventana
grid % rejilla

2. Ejercicio

Considerar la siguiente función de transferencia:

 Se pretende realizar el análisis de dicho sistema mediante el método de repuesta en frecuencia.
3. Dibujar el diagrama de Bode correspondiente para los siguientes valores de K:
 0 ≤ K < 64
 K = 64
 64 < K < 100
 K = 100
 100 < K < 260
 K = 260
 K > 260
4. Calcular los márgenes de fase y los márgenes de ganancia. Analizar la estabilidad.

5. Ejercicio

Analizar la respuesta en frecuencia de los siguientes sistemas de control con realimentación unitaria.
Dibujar los diagramas de Nyquist y de Nichols y calcular el Margen de Fase, el Margen de ganancia,
ancho de banda, pico de resonancia y frecuencia de resonancia.

V OBSERVACIONES Y CONCLUSIONES

V.1. Haga sus observaciones y emita al menos cinco conclusiones del trabajo realizado
___________________________________________________________________________________
___________________________________________________________________________________
___________________________________________________________________________________
__________________________________________________________________________________
 SOLUCIÓN LABORATORIO 5

1. Para fines didácticos, iremos analizar la respuesta en frecuencia de un sistema de suspensión de

automóvil.

%//////////////////////////////////////////////
% Sistema de Suspensión //
% /////////////////////////////////////////////
clear % borra las variables antiguas
% /////////////////////////////////////////////
% Parámetros
m=250; % masa suspensa (kg)
k=10000; % rigidez del resorte (N/m)
b=316.23; % amortecimiento (Ns/m)
% Matrices del sistema
A=[0 1;-k/m -b/m];
B=[-1;b/m];
C=[-k/m -b/m];
D=[b/m];
% Sistema en el espacio de estados
G=ss(A,B,C,D);
% Respuesta en frecuencia
% Diagrama de Bode del sistema (observe las escalas de los graficos):
bode(G)
% La escala horizontal es logaritmica, y muestra la frecuencia en rad/s.
% La escala vertical del grafico superior es lineal, y muestra la ganancia en db.
% La escala vertical del grafico inferior es lineal y muestra la fase en grados.
% Matlab asigna valores por defecto a nuestro grafico de Bode, el Matlab
% mostro que la amplitud de la ganancia esta entre -0db e +40db, la
% fase está entre -90 e 90 grados, y que la frecuencia esta entre 10 y
% 1000 rad/s. Vamos construir un vector w donde el logaritmo de la
% distancia entre un valor y otro es constante, comenzando en 10
% elevado a (-3), terminando en 10 elevado a (4), con 2000 puntos, y
% asumiendo que la unidad es rad/s:
w=logspace(-3,4,2000);
% Determinamos los vetores de magnitud y fase,
[mag,phase]=bode(G,w);
% mag y fase son arrays de 3 dimensiones, esto quiere decir que la
% función bode trabaja con sistemas multivariable (MIMO)
% Entonces, para nuestro caso, debemos solo capturar el vector correspondiente
mag=mag(1,:);
phase=phase(1,:);
% Ploteando el Bode de ganancia
subplot(2,1,1)
semilogx(w,20*log10(mag));% semilogx realiza un plot semilogaritmico
title('Diagrama de Bode');xlabel('Frecuencia (w)');ylabel('Ganancia (dB)');
axis([10^-3 10^4 -70 60]) % definimos los márgenes de la ventana
 grid % rejilla
% Ploteando el Bode de fase
subplot(2,1,2)
semilogx(w,phase);% semilogx realiza un plot semilogaritmico
xlabel('Frecuencia (w)');ylabel('Fase (°)');
axis([10^-3 10^4 -100 100]) % definimos los márgenes de la ventana
grid % rejilla

2. Ejercicio
Considerar la siguiente función de transferencia:

Se pretende realizar el análisis de dicho sistema mediante el método de repuesta en frecuencia.

 Programa plasmando la función de transferencia anterior

K1=32; %valores de K
den1=[1 0];
den2=[1 4];
denx=conv(den1,den2);
den3=[1 4 20];
den=conv(denx,den3); %Denominador
G=tf(K1,den); %Funciòn de transferencia
3. Dibujar el diagrama de Bode correspondiente para los siguientes valores de K:

o 0 ≤ K < 64
o K = 64
o 64 < K < 100
o K = 100
o 100 < K < 260
o K = 260
o K > 260

4. Calcular los márgenes de fase y los márgenes de ganancia. Analizar la estabilidad.

 Margen de ganancia
Distancia en decibelios medido desde el gráfico de bode amplitud hasta 0dB, marcado por la
frecuencia en la que la fase llega a -180º.
Positivo cuando la distancia se encuentra por debajo de los 0dB y negativa si se encuentra por
encima de los 0dB.

 Margen de Fase
Distancia en grados medido desde el gráfico de bode de fase hasta -180º, en la frecuencia en la
que la ganancia llega a 0dB.
Positivo cuando la distancia se encuentra por encima de los -180º y negativo cuando se encuentra
por debajo de los -180º.

 Para analizar el punto 3. Y 4. se terminó el programa anterior, aquí nos detalla los márgenes de
ganancia y de fase a su respectiva frecuencia en el gráfico de bode con el comando ‘margin’. Se
eligió los valores medios para K donde existe un intervalo.

K1=360; %valores de K
den1=[1 0];
den2=[1 4];
denx=conv(den1,den2);
den3=[1 4 20];
den=conv(denx,den3); %Denominador
G=tf(K1,den); %Funciòn de transferencia
margin(G) %Graficamos el bode con los màrgenes
grid

 0 ≤ K < 64 ; K=32;
 K = 64

 64 < K < 100 ; K=82

 K = 100

 100 < K < 260 ; K=180;

 K = 260

 K > 260 ; K=360;

 Análisis de los gráficos de Bode
Se puedo observar que al aumentar la ganancia con el valor K, las frecuencias de corte en 0dB y -
180º se van acercando cada vez más y según la FT el punto crítico donde ambos son 0 es en K=260,
luego de sobrepasar dicho número ambos se vuelven negativos.

5. Ejercicio
Analizar la respuesta en frecuencia de los siguientes sistemas de control con realimentación unitaria.
Dibujar los diagramas de Nyquist y de Nichols y calcular el Margen de Fase, el Margen de ganancia,
ancho de banda, pico de resonancia y frecuencia de resonancia.

 Matlab grafica el diagrama de Nyquist y Nichols siempre y cuando se defina el sistema realimentado
unitariamente; por ello utilizamos ‘feedback’ para realimentar el sistema y Matlab extraiga la
ecuación característica.

 Programa principal para el gráfico de Nyquist y Nichols según G(s)

 G=tf(num,den); %Sistema de control en lazo abierto
F=feedback(G,1);%Sistema realimentado unitariamente
subplot(1,2,1);
nyquist(F),grid;
subplot(1,2,2);
nichols(F),grid;

 Para calcular los márgenes de ganancia y fase, pico y frecuencia de resonancia utilizamos un poco
la vista y las características ofrecidas por el gráfico. Se tomará este punto como ejemplo:

a)

clear
num=[0.1 1];
den1=[1 0];
den2=[1 1];
den3=[0.01 1];
denx=conv(den1,den2);
den=conv(denx,den3);
G=tf(num,den); %Sistema de control en lazo abierto
F=feedback(G,1);%Sistema realimentado unitariamente
subplot(1,2,1);
nyquist(F),grid;
subplot(1,2,2);
nichols(F),grid;

 Para visualizar el pico y frecuencia de resonancia seleccionamos ‘peak response’ en ambas

gráficas para confirmar.
 Pico de resonancia (dB) = 0.79
Frecuencia de resonancia (rad/s) = 0.61

 Para visualizar el margen de ganancia y fase seleccionamos ‘Minimun Stability Margins’ en ambas
gráficas para confirmar.
 Para Nyquist, el margen de fase mide desde el punto (-1+j0) hasta la intersección de la magnitud
1 y la gráfica; si el ángulo se mide hacia abajo, el margen de fase es positiva, si se mide hacia
arriba, es negativa; todo siempre y cuando se evalúe en frecuencia positiva.

 Y el margen de ganancia se calcula midiendo L desde (0,j0) hasta la intersección de la gráfica en

1 1
la parte izquierda del eje imaginario. Si 𝐿
> 1 el margen de ganancia es positiva, y si 𝐿
< 1 es
1
negativa. Para obtener en decibelios se transforma con 𝑑𝐵 = 20 ∗ log (𝐿) .

 Margen de ganancia y fase pueden observarse mejor en Nichols.

Margen de Ganancia (dB) = +Infinito
Margen de Fase (deg) = +105

 Y por último para el ancho de banda se hace un recorrido en la gráfica de Nichols hasta cuando la
ganancia baje a -3dB. En otras palabras, el ancho de banda se encuentra es la gamma de
frecuencias donde la ganancia es superior a -3dB.
 Ancho de Banda (rad/s) = [0 , 1.22]

b)
clear
num1=[0.5];
num2=[1 1];
num=conv(num1,num2);
den1=[1 0];
den2=[0.2 1];
den3=[0.5 1 1];
denx=conv(den1,den2);
den=conv(denx,den3);
G=tf(num,den); %Sistema de control en lazo abierto
F=feedback(G,1);%Sistema realimentado unitariamente
subplot(1,2,1);
nyquist(F),grid;
subplot(1,2,2);
nichols(F),grid;
 Pico de Resonancia (dB) = 0
Frecuencia de Resonancia (rad/s) = 0
Margen de Ganancia (dB) = +15.8
Margen de Fase (deg) = -180
Ancho de Banda (rad/s) = [0 , 0.993]

c)

clear
num=[1 1];
den1=[1 0];
den2=[0.2 1];
den3=[0.5 1];
denx=conv(den1,den2);
den=conv(denx,den3);
G=tf(num,den); %Sistema de control en lazo abierto
F=feedback(G,1);%Sistema realimentado unitariamente
subplot(1,2,1);
nyquist(F),grid;
subplot(1,2,2);
nichols(F),grid;
 Pico de Resonancia (dB) = 0
Frecuencia de Resonancia (rad/s) = 0
Margen de Ganancia (dB) = +infinito
Margen de Fase (deg) = -180
Ancho de Banda (rad/s) = [0 , 1.16]

d)
clear
num=[1];
den1=[1 0];
den2=[1 1];
den3=[0.5 1];
denx=conv(den1,den2);
den=conv(denx,den3);
G=tf(num,den); %Sistema de control en lazo abierto
F=feedback(G,1);%Sistema realimentado unitariamente
subplot(1,2,1);
nyquist(F),grid;
subplot(1,2,2);
nichols(F),grid;
 Pico de Resonancia (dB) = 5.28
Frecuencia de Resonancia (rad/s) = 0.824
Margen de Ganancia (dB) = +6.03
Margen de Fase (deg) = +24.2
Ancho de Banda (rad/s) = [0 , 1.26]

e)
clear
num=[50];
den1=[1 0];
den2=[1 1];
den3=[0.5 0 1];
denx=conv(den1,den2);
den=conv(denx,den3);
G=tf(num,den); %Sistema de control en lazo abierto
F=feedback(G,1);%Sistema realimentado unitariamente
subplot(1,2,1);
nyquist(F),grid;
subplot(1,2,2);
nichols(F),grid;
 Pico de Resonancia (dB) = 0.0864
Frecuencia de Resonancia (rad/s) = 1.04
Margen de Ganancia (dB) = infinito
Margen de Fase (deg) = +180
Ancho de Banda (rad/s) = [0 , 1.26]

f)

clear
num1=[0.2];
num2=[-1 0 1];
num3=[0.5 1];
numx=conv(num1,num2);
num=conv(numx,num3);
den1=[1 0];
den2=[1 1 1];
den=conv(den1,den2);
G=tf(num,den); %Sistema de control en lazo abierto
F=feedback(G,1);%Sistema realimentado unitariamente
subplot(1,2,1);
nyquist(F),grid;
subplot(1,2,2);nichols(F),grid;
 Pico de Resonancia (dB) = 0
Frecuencia de Resonancia (rad/s) = 0
Margen de Ganancia1 (dB) = +7.37
Margen de Ganancia2 (dB) = +19.1
Margen de Fase (deg) = -180
Ancho de Banda (rad/s) = [0 , 0.256] , [0.891 , 0.997]

VII OBSERVACIONES Y CONCLUSIONES

 Se ha podido observar los márgenes de ganancia y fase para funciones de transferencia y sistemas
realimentados.
 Se ha entendido como trabaja Matlab los diagramas de Bode, Nyquist y Nichols, los cuales se ha
comentado anteriormente en la realización de este laboratorio.
 Se ha aplicado un poco de lógica para obtener ciertas magnitudes como lo es el margen de
ganancia viendo la gráfica.
 Comparando, el diagrama para observar mejor el comportamiento en frecuencia es Bode.
 El diagrama de Nichols es muy similar al de bode lo cual hace que obtener las características de
la gráfica no sea tan complicado.
 El diagrama de Nyquist que se grafica en los ejes real e imaginario es un poco complicado de
observar ya que generalmente se trabaja en grados y decibelios para analizar funciones de
transferencia, pero de igual manera se logra apreciar las características interpretando y
convirtiendo.
 Hemos podido comprender un poco sobre la estabilidad de los sistemas realimentados, por
ejemplo: del inciso 5, de a-e, concluimos que un sistema realimentado es estable siempre y cuando
sus márgenes de ganancia y fase sean positivos. El punto f del inciso 5, es un poco mas complicado
de comprender debido a que su numerador es de un orden superior a los trabajados anteriormente
lo cual nos produjo doble margen de ganancia y doble ancho de banda.