• Author / Uploaded
• Kevin Medina Villagomez

• Categories
• Verdad lógica
• Teoría de sistemas
• Algoritmos
• Física y matemáticas
• Matemáticas

Citation preview

UNIVERSIDAD NACIONAL DE SAN AGUSTIN Facultad de Ingeniería de Producción y Servicios Escuela Profesional de Ingeniería Electrónica

Control Automatico 1 LABORATORIO N° 2 “SIMULACION NUMERICA”

Arequipa - Perú 2015

PRÁCTICA NÚMERO 2 SIMULACION NUMERICA 1. Ecuaciones diferenciales (repetir la siguiente actividad). Determine la ecuación diferencial y’=[1-y(t)]/2 usando el método de Euler. Plotear el grafico de y(t), en función de t, y el grafico de la solución exacta ye(t) en función de t. Instante inicial: t=0 Instante final: tf=10 Solución exacta de la ecuación diferencial: ye(t) =1-e^(1-t/2) function[ydot]=edo1(y) ydot=(1-y)/2

%definir una funcion

t(1)=0; tf=10; y(1)=0; ye(1)=0; h=0.5; n=round(tf/h); for i=1:n t(i+1)=t(i)+h; ydot(i)=edo1(y(i)); y(i+1)=y(i)+h*ydot(i); ye(i+1)=1-exp(-t(i+1)/2); end ydot(n+1)=edo1(y(n+1)); plot(t,y,t,ye,'r:'); legend('solución exacta','solución numerica-euler') title('comparacion entre la solución exacta y la solución numerica') xlabel('tiempo(s)');ylabel('amplitud()') euler_ode1 comparacion entre la solucion exacta y la solucion numerica

1

solucion exacta solucion numerica-euler

0.9 0.8 0.7 amplitud()

 

0.6 0.5 0.4 0.3 0.2 0.1 0

0

1

2

3

4

5 tiempo(s)

6

7

8

9

10

2. Repita las instrucciones: F.T. H(s)=1/(s+1) h=tf(1,[1 1])

g=ss([0 1;-5 -2],[0;3],[0 1],0)

[num,den]=ss2tf([0 1;-5 -2],[0;3],[0 1],0)

g=tf(16,[1 9 16]) subplot(2,2,1) bode(g) subplot(2,2,2) pzmap(g) subplot(2,2,3) step(g) subplot(224)

Bode Diagram

0 -50 -100 0 -90 -180 -1 10

0

1

10

Pole-Zero Map

1 Imaginary Axis

Phase (deg)Magnitude (dB)

nyquist(g)

0 -0.5 -1 -8

2

10

0.5

10

-6

-4

Step Response

1 Imaginary Axis

Amplitude

1

0.5

0

0

0.5

1

1.5

2

2.5

-2

0

Real Axis

Frequency (rad/sec)

Nyquist Diagram

0.5 0 -0.5 -1 -1

-0.5

Time (sec)

Diagrama: de bodes, Nyquis

0 Real Axis

0.5

1

3. Ejercicio: Realice la simulación numérica de la ecuación diferencial descrita en el item1; Esta vez implemente el algoritmo RK4, en ves del método euler. Haga una comparación entre las soluciones: exacta, euler y RK4. function[f]=rk(x,y); f=(1-y)/2;

x(1)=0; tf=10 y(1)=0; ye(1)=0; h=0.5; n=round(tf/h); for

i=1:n;

k1=f(x(i),y(i)); k2=f(x(i)+h/2,y(i)+h/2*k1); k3=f(x(i)+h/2,y(i)+h/2*k2); k4=f(x(i)+h,y(i)+h*k3); y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4); x(i+1)=x(i)+h; ye(i+1)=1-exp(-t(i+1)/2); end ydot(n+1)=rk(y(n+1)); plot(x,y,x,ye,'r:'); legend('solución exacta','solución numerica-Euler'); title('Comparacion entre la solución exacta y la solución numerica') xlabel('Tiempo(s)') ylabel('Amplitud( )')

4. Ejercicio: Repetir uno , 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 :

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

t(1)=0; tf=600; y(1)=0.9; ye(1)=0.9; h=0.85; n=round(tf/h);

% Funcion de transferencia

for i=1:n t(i+1)=t(i)+h; ydot(i)=edo1(y(i)); y(i+1)=y(i)+h*ydot(i); ye(i+1)=1-exp(-t(i+1)/2); end ydot(n+1)=edo1(y(n+1)); plot(t,y,t,ye,'r:'); legend('solución exacta','solución numerica-euler') title('comparacion entre la solución exacta y la solución numerica') xlabel('tiempo(s)');ylabel('amplitud()')

euler_ode1 comparacion entre la solucion exacta y la solucion numerica

1.3

solucion exacta solucion numerica-euler

1.2 1.1

amplitud()

1 0.9 0.8 0.7 0.6 0.5 0.4 0

100

200

300 400 tiempo(s)

500

600

700

5. Ejercicio: reducir el diagrama de bloques de la figura 2 y expresar la F.T. entre C y R en función de ( G1,G2 y H ). G1

+ R

+

G2

C

+

H Solución:

G1 +

R

+

G2

+

+ C

H G1 G2

R

+

+

+

G2

C

+ H

G2 1+G2*H

R

R

1+G1 G2

G1+G2 1+G2*H

C

C

6. Ejercicio: el propósito de estos ejercicios es familiarizarse con las funciones de transferencia, su representación y respuesta. Ud. Podrá explorar la relación entre la localización de polos y zeros, la respuesta temporal y la respuesta frecuencial.Por tanto, para las siguientes funciones de transferencia grafique (i)el diagrama de polos y zeros (ii)el diagrama de bode y (iii)la respuesta escalón unitario. a) Sistema de primer orden: H1(s)=

1 s 1

Pole-Zero Map

Imaginary Axis

1 0.5 0 -0.5 -1 -1

-0.5

0

Real Axis

Amplitude

0.5

0

0

2 Time (sec)

0

Bode Diagram

-20 -40 0 -45 -90 -2 10

0

10

Frequency (rad/sec)

Step Response

1

Phase (deg) Magnitude (dB)

g=tf(1,[1 1]); % función de transferencia subplot(2,2,1); % divide la figura en 2 filas, 2 columnas y ubica la grafica en la posición 1 pzmap(g); % grafica de los polos y ceros de la función de transferencia subplot(2,2,2); % divide la figura y ubica la grafica en la posición 2 bode(g); % bode de amplitud y fase de la función de transferencia subplot(2,2,3); % divide la figura y ubica la grafica en la posición 3 step(g); % gráfica de la respuesta escalón

4

6

10

2

b) Sistema de segundo orden:

H2(s)=

2  s  1 s  2

Pole-Zero Map

Imaginary Axis

1 0.5 0 -0.5 -1 -2

-1.5

-1

-0.5

0

Real Axis

Amplitude

0.5

0

0

5 Time (sec)

0

Bode Diagram

-50 -100 0 -90 -180 -1 10

0

10

1

10

Frequency (rad/sec)

Step Response

1

Phase (deg) Magnitude (dB)

g=tf(2,[conv([1 1],[1 2])]); % función de transferencia subplot(2,2,1); % divide la figura en 2 filas, 2 columnas y ubica la grafica en la posición 1 pzmap(g); % grafica de los polos y ceros de la función de transferencia subplot(2,2,2); % divide la figura y ubica la grafica en la posición 2 bode(g); % bode de amplitud y fase de la función de transferencia subplot(2,2,3); % divide la figura y ubica la grafica en la posición 3 step(g); % gráfica de la respuesta escalón

10

2

10

c) Sistema de primer orden (adicionamos un zero, primero en el RHP, luego en el LHP):

H3(s)=

1 s2 2  s  1

Pole-Zero Map

Imaginary Axis

1 0.5 0 -0.5 -1 -1

0

1

2

Real Axis

Amplitude

0 -0.5 -1

0

2 Time (sec)

0

Bode Diagram

-5 -10 180 90 0 -1 10

0

10

1

10

Frequency (rad/sec)

Step Response

0.5

Phase (deg)Magnitude (dB)

g=tf(0.5*[1 -2],[1 1]); % función de transferencia subplot(2,2,1); % divide la figura en 2 filas, 2 columnas y ubica la grafica en la posición 1 pzmap(g); % grafica de los polos y ceros de la función de transferencia subplot(2,2,2); % divide la figura y ubica la grafica en la posición 2 bode(g); % bode de amplitud y fase de la función de transferencia subplot(2,2,3); % divide la figura y ubica la grafica en la posición 3 step(g); % gráfica de la respuesta escalón

4

6

2

10

H4(s)=

1 s2 2  s  1

Pole-Zero Map

Imaginary Axis

1 0.5 0 -0.5 -1 -2

-1.5

-1

-0.5

0

Real Axis

Amplitude

0.8 0.6 0.4

0

2 Time (sec)

0

Bode Diagram

-5 -10 0 -10 -20 -1 10

0

10

1

10

Frequency (rad/sec)

Step Response

1

Phase (deg)Magnitude (dB)

g=tf(0.5*[1 2],[1 1]); % función de transferencia subplot(2,2,1); % divide la figura en 2 filas, 2 columnas y ubica la grafica en la posición 1 pzmap(g); % grafica de los polos y ceros de la función de transferencia subplot(2,2,2); % divide la figura y ubica la grafica en la posición 2 bode(g); % bode de amplitud y fase de la función de transferencia subplot(2,2,3); % divide la figura y ubica la grafica en la posición 3 step(g); % gráfica de la respuesta escalón

4

6

2

10

d) El siguiente sistema es un sistema de segundo orden, con Wn=5 y un ζ=0.2 25 H(5)= -------------s^2 + 2 s + 25

Pole-Zero Map

Imaginary Axis

5

0

-5 -1

-0.5

0

Real Axis

Amplitude

1.5 1 0.5 0

0

2

4 Time (sec)

100

Bode Diagram

0 -100 0 -90 -180 -1 10

0

10

1

10

Frequency (rad/sec)

Step Response

2

Phase (deg)Magnitude (dB)

g=tf(25,[1 2 25]); % función de transferencia subplot(2,2,1); % divide la figura en 2 filas, 2 columnas y ubica la grafica en la posición 1 pzmap(g); % grafica de los polos y ceros de la función de transferencia subplot(2,2,2); % divide la figura y ubica la grafica en la posición 2 bode(g); % bode de amplitud y fase de la función de transferencia subplot(2,2,3); % divide la figura y ubica la grafica en la posición 3 step(g); % gráfica de la respuesta escalón

6

2

10

e) Ahora adicionamos un cero sobre el origen: 25 s H6 = ----------------s^2 + 2 s + 25

Pole-Zero Map

Imaginary Axis

5

0

-5 -1

-0.5

0

Real Axis

Amplitude

2 0 -2

0

2 Time (sec)

50

Bode Diagram

0 -50 90 0 -90 -1 10

0

10

1

10

Frequency (rad/sec)

Step Response

4

Phase (deg)Magnitude (dB)

g=tf([25 0],[1 2 25]); % función de transferencia subplot(2,2,1); % divide la figura en 2 filas, 2 columnas y ubica la grafica en la posición 1 pzmap(g); % grafica de los polos y ceros de la función de transferencia subplot(2,2,2); % divide la figura y ubica la grafica en la posición 2 bode(g); % bode de amplitud y fase de la función de transferencia subplot(2,2,3) % divide la figura y ubica la grafica en la posición 3 step(g); % gráfica de la respuesta escalón

4

6

2

10

e) Ahora adicionamos un cero a 5 rad/seg:

H7=

5 s + 25 -------------s^2 + 2 s + 25

g=tf(5*[1 5],[1 2 25]); % función de transferencia subplot(2,2,1); % divide la figura en 2 filas, 2 columnas y ubica la grafica en la posición 1 pzmap(g); % grafica de los polos y ceros de la función de transferencia subplot(2,2,2); % divide la figura y ubica la grafica en la posición 2 bode(g); % bode de amplitud y fase de la función de transferencia subplot(2,2,3); % divide la figura y ubica la grafica en la posición 3 step(g); % gráfica de la respuesta escalón

0

-5 -6

-4

-2

0

Real Axis

Amplitude

1.5 1 0.5 0

0

2 Time (sec)

50

Bode Diagram

0 -50 180 0

-180 -1 10

0

10

1

10

Frequency (rad/sec)

Step Response

2

Phase (deg)

Imaginary Axis

Magnitude (dB)

Pole-Zero Map

5

4

6

2

10

f) Un tercer polo es adicionado a H(5), primero a baja frecuencia, y luego a alta frecuencia; 25 H8 = ---------------------------s^3 + 3 s^2 + 27 s + 25

y

H9 =

125 ---------------------------s^3 + 7 s^2 + 35 s + 125

Bode Diagram

0

5

-50

Imaginary Axis

Phase (deg) Magnitude (dB)

g1=tf(25,conv([1 1],[1 2 25])) ; % función de transferencia 1 g2=tf(125,conv([1 5],[1 2 25])); % función de transferencia 2 figure(1); % figura numero 1 subplot(2,2,1); % divide la figura 2 filas y columnas ubicado en 1 bode(g1); % bode de amplitud y fase de la función de transferencia subplot(2,2,2); % divide la figura y la ubica en la posición 2 pzmap(g1); % grafica de polos y ceros de la función de transferencia subplot(2,2,3); % divide la figura y la ubica en la posición 3 step(g1); % grafica de la respuesta escalón subplot(2,2,4); % divide la figura y la ubica en la posición 4 nyquist(g1); % grafica de Nyquist de la función de transferencia figure(2); % figura numero 2 subplot(2,2,1); % divide la figura en 2 filas y columnas ubicado en 1 bode(g2); % bode de amplitud y fase de la función de transferencia subplot(2,2,2); % divide la figura y la ubica en la posición 2 pzmap(g2); % grafica los polos y ceros de la función de transferencia subplot(2,2,3); % divide la figura y la ubica en la posición 3 step(g2); % grafica de la respuesta escalón subplot(2,2,4); % divide la figura y la ubica en la posición 4 nyquist(g2); % grafica de Nyquist de la función de transferencia

-100 0 -180 -360 -2 10

0

0

-5 -1.5

2

10

10

Pole-Zero Map

-1

Frequency (rad/sec) Step Response

1 Imaginary Axis

Amplitude

1

0.5

0

0

2 Time (sec)

4

6

-0.5

0

Real Axis Nyquist Diagram

0.5 0 -0.5 -1 -1

-0.5

0 Real Axis

0.5

1

5

0

Imaginary Axis

Phase (deg)Magnitude (dB)

Bode Diagram

100

-100 0 -180 -360 -1 10

0

1

10

0

-5 -6

2

10

10

Pole-Zero Map

-4

Frequency (rad/sec) Step Response

2 Imaginary Axis

Amplitude

1.5 1 0.5 0

0

1

2

3

Time (sec)

4

5

-2

0

Real Axis Nyquist Diagram

1 0 -1 -2 -2

-1

0

1

Real Axis

CONCLUSIONES  Con el metodo de Runge Kutta, optimizamos el metodo de Euler, para obtener soluciones mucho más exactas  Usando los comandos que nos proporcional el matlab podemos obtener con mayor facilidad los diagramas de bode, Nyquist, la respuesta escalón, pero para ello debemos saber reducir cualquier diagrama de bloques o grafos para poder reducir todo a una sola función de transferencia, para ello observamos las distintas formas que hay en el libro de ogata.