MODELACION DE PROCESOS

5/10/07

Problema Nº1

La Fig. Nº1, muestra un levitador magnético idealizado. Donde la tensión e(t) es manipulada para que la
bolita se sitúe a 30 cm sobre el nivel del suelo. Para esto se utiliza el electroimán que puede ser representado
por el circuito R L. Recordar que la dinámica impuesta por este circuito, está dada por una ecuación
diferencial de 1º orden, proveniente de LVK Ley de Voltaje de Kirchoff, desarrollada en (1).

(1)

Por otro lado, la fuerza de atracción Fe(t), que ejerce el electroimán, está dada por (2).

, (2)

Unidad de
Parámetro Valor
Medida
M 0.2

d 1.5

k 24.5

g 9.8

l0 0,3
l1 0.6

ki 3e−3

a 0.02
R 1

L 0.06

1
Fig. Nº1 Levitador magnético.

a. Determinar las variables de entrada, variables de salida, variables de estado, perturbaciones y parámetros.

b. Escribir las ecuaciones asociadas al sistema metódicamente.

c. Determine la tensión e10, e20, y e30 a aplicar en t0 = 0, t1 = 1, t2 = 5 s, para tener en S.S. (Stationary State)
a la bolita a 30, 35 y 25 cm desde el piso, respectivamente.

d. Encuentre el perfil de tensión a aplicar para conseguir lo indicado en anteriormente y simule el sistema para
un rango de simulación de 0 " t < 8 s. Grafique e[V], i[A], x[cm], y v [km/hr]. Asegúrese de estar en S.S. para
t = 0.

e. Comentar los resultados obtenidos.

Solución Tarea N°1 Modelación de Procesos.

Problema N°1

a. Determinar las variables de entrada, variables de salida, variables de estado, perturbaciones y

parámetros.

Variables de entrada:

• Tensión: e(t)

Variables de salida:

• Posición: x(t)

Variables de estado:

• Corriente del electroimán: i(t)

• Posición: x(t)
• Velocidad:

2
Perturbaciones:

• No hay

Parámetros:

• Masa: M
• Coeficiente de viscosidad del amortiguador: d
• Coeficiente de elasticidad del resorte: k
• Aceleración de gravedad: g
• Distancia de la bolita al piso: l0
• Distancia del electroimán al piso: l1
• Constante de corriente del electroimán: ki
• Espesor del imán: a
• Resistencia: R
• Inductancia de la bobina: L

b. Escribir las ecuaciones asociadas al sistema metódicamente.

Para este modelo se han tomado ciertas consideraciones para que sea un levitador magnético idealizado, por
ejemplo:

• No existe fricción viscosa entre la masa levitante y el aire.

• La masa M a levitar solo tiene un grado de libertad (x).

A partir de la posición de equilibrio la fuerza de equilibrio, la fuerza que realiza el resorte puede escribirse
como:

La fuerza que realiza el amortiguador como:

La fuerza de gravedad:

La fuerza del electroimán:

Figura 2. Diagrama de cuerpo libre sistema

3
Masa−Resorte−Amortiguador−electroimán

De la segunda Ley de Newton y diagrama de cuerpo libre de la Figura 2, se puede escribir:

Es decir:

, (3)

Tenemos que consideras que:

, (4)

Y así dejamos y(t) en función de x(t).

Reemplazando la ecuación (4) en (3) y ordenando, nos queda:

(5)

También podemos obtener otra ecuación de la ecuación (1):

(1)

Despejando

, tenemos la ecuación (6):

(6)

Ocuparemos (1) y (5) para representar nuestro sistema en matrices, con

,

4
,

y
:

Claramente se ve que el sistema de ecuaciones es no lineal. Para linealizarlo encontramos sus derivadas con
respecto a las variables de estado, entrada y perturbaciones para obtener el modelo.

Ahora tenemos el siguiente sistema de ecuaciones LINEAL:

c. Determine la tensión e10, e20, y e30 a aplicar en t0 = 0, t1 = 1, t2 = 5 s, para tener en S.S. (Stationary
State) a la bolita a 30, 35 y 25 cm desde el piso, respectivamente.

Sabemos por las ecuación (1) que e depende de i, por lo tanto basta encontrar i en el tiempo indicado para
saber que voltaje aplicar. Ocupamos la ecuación (5), sabiendo que para que se cumpla la condición de sistema
estacionario, no debe haber aceleración y la velocidad tiene que ser una constante, así:

5
Ahora solo queda reemplazar la posición y los parámetros en la ecuación, dando i(t)=cte, reemplazándola en
la ecuación (1), teniendo los siguientes valores:

e10= 14.46 [V]

e20= 16.95 [V]

e30= 9.52 [V]

d. Encuentre el perfil de tensión a aplicar para conseguir lo indicado en anteriormente y simule el

sistema para un rango de simulación de 0 " t < 8 s. Grafique e[V], i[A], x[cm], y v [km/hr]. Asegúrese de
estar en S.S. para t = 0.

Los graficos se obtuvieron gracias a los siguientes códigos:

−E1.m:

%Tarea 1: PROBLEMA N°1: SISTEMA MASA RESORTE AMORTIGUADOR ELECTROIMAN

clear all;

close all;

clc;

global M d k g l0 l1 ki a R L

%Parámetros

M=0.2; %[Kg]

d=1.5; %[Kg/s]

k=24.5; %[Kg/s^2]

g=9.8; %[m/s^2]

6
l0=0.3; %[m]

l1=0.6; %[m]

ki=3e−3;%[Kg*m^2/(A^2*s^2)]

a=0.02; %[m]

R=1; %[ohm]

L=0.06; %[H]

%Condiciones iniciales:

i0=14.46;%[A]

x0=0.3; %[m]

v0=0; %[m/s]

%Problemas Stiff

[t,x]=ode23s('E1_ec',[0 8],[i0 x0 v0]');

U=14.46+2.47*(t>1)−9.52*(t>5);

figure(1)

subplot(311),plot(t,x(:,1),'r')

grid on, xlabel('tiempo (s)'),ylabel('Corriente [A]'),title('i [A]')

subplot(312),plot(t,x(:,2)*1e2,'g')

grid on, xlabel('tiempo (s)'), ylabel('Posición [cm]'),title('x [cm]')

subplot(313),plot(t,x(:,3)*3.6)

grid on, xlabel('tiempo (s)'), ylabel('Velocidad [Km/hr]'),title('v [Km/hr]')

figure(2)

plot(t,U,'r'), axis([0 8 0 20])

grid on, xlabel('tiempo (s)'), ylabel('Tensión [V]'),title('e [V]')

y E1_ec.m:

function v_estado=E1_ec(t,x)

global M d k g l0 l1 ki a R L

7
U=14.46+2.47*(t>1)−9.52*(t>5);

v_estado=zeros(3,1);

v_estado(1) = −(R/L)*x(1) + U/L;

v_estado(2) = x(3);

v_estado(3) = −g − k*(x(2)−l0)/M − d*x(3)/M + ki*(x(1)^2)/((l1−x(2)+a)*M);

e. Comentar los resultados obtenidos.

Para obtener estado estacionario en t=0 se simulo un intervalo de 0 a 8 segundos. De este modo podemos
apreciar el movimiento de la bolita, que en cada cambio de tensiones tiene su estado estacionario en tiempos
cercanos a 1 segundo, lo que es bastante rápido, esto gracias a la acción del amortiguador. Podemos notar que
el tiempo que tarda en obtener estado estacionario es proporcional a la diferencia de tensión que se le a
aplicado.

Problema Nº2

El brazo robótico SCARA, mostrado en la Fig. Nº2, es un importante manipulador, el cual posee 4 grados de
libertad (DOF, Degree Of Freedom), sin considerar el gripper. Es producido por varias compañías, entre las
cuales se puede nombrar Adept Technology, IBM, Seiko, entre otras. El problema a desarrollar, consiste en
simular en 3D, mediante Matlab®, este manipulador robótico, donde las variables de entrada al programa,
serán los ángulos , , y la altura lineal d3. Como ayuda adicional son incorporadas las matrices de
rotación y traslación según corresponda. Es tarea del lector verificar si están bien escritas.

Unidad de
Símbolo Valor
Medida
d1 0.5
d4 1.5
a1 0.4
a2 0.35

8
Problema N°2

Variables de entrada:

• Ángulo de giro Link 1 con respecto a Joint 1:

• Ángulo de giro Link 2 con respecto a Joint 2:
• Ángulo de giro de Link 4:
• Distancia de x2 a x3: d3

Parámetros:

• d1=50 [cm]
• d4=15 [cm]

9
 • Largo de Link 1 a1=40 [cm]
• Largo de Link 2 a2=35 [cm]

Código MATLAB:

clear all;

close all;

clc;

%Entradas

t1 =input('Ingrese Ángulo Base (Link 1) (theta1[grad]) =');

t2 =input('Ingrese Ángulo Hombro(Link 2) (theta2[grad]) =');

t4 =input('Ingrese Ángulo Codo2 (Link 4) (theta4[grad]) =');

d3 =input('Ingrese distancia de x2 a x3 (d3[cm]) =') ;

%==========================================================================

%==========================================================================

% Traspasa ángulo a radianes:

theta1 =t1*pi/180;

theta2 =t2*pi/180;

theta4 =t4*pi/180;

%==========================================================================

%==========================================================================

%Parámetros del modelo:

d1=50; %Distancia de la base al entro de Link 1 [cm]

d4=15; %Largo de Link 4 [cm]

a1=40; %Largo de Link 1 [cm]

a2=35; %Largo de Link 2 [cm]

%==========================================================================

%==========================================================================

% Matriz de Link 1 respecto a Link 0 [ESTABA BUENA]

10
A1 = [ cos(theta1) −sin(theta1) 0 a1*cos(theta1);

sin(theta1) cos(theta1) 0 a1*sin(theta1);

0 0 1 d1 ;

0 0 0 1] ;

%==========================================================================

%==========================================================================

% Matriz de Link 2 respecto a Link 1 [ESTABA MALA]

A2 = [ cos(theta2) sin(theta2) 0 a2*cos(theta2);

sin(theta2) −cos(theta2) 0 a2*sin(theta2);

0010;

0 0 0 1] ;

%OBS:

% Matriz que había en el enunciado original:

% A2 =[cos(theta2) sin(theta2) 0 a2*cos(theta2);

% sin(theta2) −cos(theta2) 0 a2*sin(theta2);

% 0 0 "−1" 0 ; −1 reemplazado por 1

% 0 0 0 1] ;

%==========================================================================

%==========================================================================

% Matriz que baja la posición [ESTABA MALA]

A3 = [ 1 0 0 0 ;

0100;

0 0 1 −d3;

0 0 0 1];

%OBS:

11
% Matriz que había en el enunciado original:

% A3 = [ 1 0 0 0 ;

%0100;

% 0 0 1 "d3"; d3 reemplazado por −d3

% 0 0 0 1] ;

%==========================================================================

%==========================================================================

%Matriz giro de Link 4

A4 = [ cos(theta4) −sin(theta4) 0 0 ;

sin(theta4) cos(theta4) 0 0 ;

0 0 1 −d4;

0 0 0 1];

% Matriz que había en el enunciado original:

% A4 = [ cos(theta4) −sin(theta4) 0 0 ;

% sin(theta4) cos(theta4) 0 0 ;

% 0 0 1 "d4" ; d4 reemplazado por −d4

% 0 0 0 1] ;

%==========================================================================

%==========================================================================

% Cálculo de vector de posición.

P0= [1 0 0 0;

0 1 0 0;

0 0 1 0;

0 0 0 1];

PR= [1 0 0 0;

0 1 0 0;

12
0 0 1 d1;

0 0 0 1];

P1= A1;

P2= A1*A2;

P3= A1*A2*A3;

P4= A1*A2*A3*A4;

%==========================================================================

%==========================================================================

%Gráfico:

x=[ P0(1,4) PR(1,4) P1(1,4) P2(1,4) P3(1,4) P4(1,4) ];

y=[ P0(2,4) PR(2,4) P1(2,4) P2(2,4) P3(2,4) P4(2,4) ];

z=[ P0(3,4) PR(3,4) P1(3,4) P2(3,4) P3(3,4) P4(3,4) ];

axis on

grid on

axis([ −100 100 −100 100 0 100]);

view([112.5,30]), hold

plot3(x,y,z)

xlabel('Eje X')

ylabel('Eje Y')

zlabel('Eje Z')

% plot del eje x0,y0,z0

x0=P0*[ 0.1*a1; 0; 0; 1];

y0=P0*[ 0; 0.1*a1; 0; 1];

z0=P0*[ 0; 0; 0.1*a1; 1];

plot3([P0(1,4) x0(1)],[P0(2,4) x0(2)],[P0(3,4) x0(3)],'r');

plot3([P0(1,4) y0(1)],[P0(2,4) y0(2)],[P0(3,4) y0(3)],'g');

13
plot3([P0(1,4) z0(1)],[P0(2,4) z0(2)],[P0(3,4) z0(3)],'k');

% plot del eje xR,yR,zR

xR=PR*[ 0.1*a1; 0; 0; 1];

yR=PR*[ 0; 0.1*a1; 0; 1];

zR=PR*[ 0; 0; 0.1*a1; 1];

plot3([PR(1,4) xR(1)],[PR(2,4) xR(2)],[PR(3,4) xR(3)],'r');

plot3([PR(1,4) yR(1)],[PR(2,4) yR(2)],[PR(3,4) yR(3)],'g');

plot3([PR(1,4) zR(1)],[PR(2,4) zR(2)],[PR(3,4) zR(3)],'k');

% plot del eje x1,y1,z1

x1=P1*[ 0.1*a1; 0; 0; 1];

y1=P1*[ 0; 0.1*a1; 0; 1];

z1=P1*[ 0; 0; 0.1*a1; 1];

plot3([P1(1,4) x1(1)],[P1(2,4) x1(2)],[P1(3,4) x1(3)],'r');

plot3([P1(1,4) y1(1)],[P1(2,4) y1(2)],[P1(3,4) y1(3)],'g');

plot3([P1(1,4) z1(1)],[P1(2,4) z1(2)],[P1(3,4) z1(3)],'k');

% plot del eje x2,y2,z2

x2=P2*[ 0.1*a2; 0; 0; 1];

y2=P2*[ 0; 0.1*a2; 0; 1];

z2=P2*[ 0; 0; 0.1*a2; 1];

plot3([P2(1,4) x2(1)],[P2(2,4) x2(2)],[P2(3,4) x2(3)],'r');

plot3([P2(1,4) y2(1)],[P2(2,4) y2(2)],[P2(3,4) y2(3)],'g');

plot3([P2(1,4) z2(1)],[P2(2,4) z2(2)],[P2(3,4) z2(3)],'k');

% plot del eje x3,y3,z3

x3=P3*[ 0.1*a2; 0; 0; 1];

y3=P3*[ 0; 0.1*a2; 0; 1];

z3=P3*[ 0; 0; 0.1*a2; 1];

14
plot3([P3(1,4) x3(1)],[P3(2,4) x3(2)],[P3(3,4) x3(3)],'r');

plot3([P3(1,4) y3(1)],[P3(2,4) y3(2)],[P3(3,4) y3(3)],'g');

plot3([P3(1,4) z3(1)],[P3(2,4) z3(2)],[P3(3,4) z3(3)],'k');

% plot del eje x4,y4,z4

x4=P4*[ 0.1*a2; 0; 0; 1];

y4=P4*[ 0; 0.1*a2; 0; 1];

z4=P4*[ 0; 0; 0.1*a2; 1];

plot3([P4(1,4) x4(1)],[P4(2,4) x4(2)],[P4(3,4) x4(3)],'r');

plot3([P4(1,4) y4(1)],[P4(2,4) y4(2)],[P4(3,4) y4(3)],'g');

plot3([P4(1,4) z4(1)],[P4(2,4) z4(2)],[P4(3,4) z4(3)],'k');

Problema Nº3

El sistema consiste en un movimiento de un proyectil en dos dimensiones, para una velocidad inicial VO, un
ángulo de elevación , además se agrega el efecto de la resistencia del aire, esta fuerza es causada por el aire
que es directamente proporcional a la velocidad del proyectil.

Se pide desarrollar detalladamente lo siguiente:

a. Determinar el modelo dinámico del sistema.

b. Obtener de manera analítica las coordenadas X e Y

c. Obtener una expresión aproximada para el tiempo de vuelo y la distancia recorrida.

d. Simule el sistema, obtenga gráficas para el desplazamiento [km] y velocidad [m/s], para el modelo obtenido
en a. y b. , con los parámetros dados en la tabla, y compare.

e. Obtenga las trayectorias para los siguientes valores de k, con k= 0.005, 0.01, 0.02, 0.04 y 0.08.

f. De la expresión aproximada en c. , para el tiempo de vuelo, cual es el rango de valor de k, para que la
aproximación sea válida.

Unidad de
Símbolo Valor
Medida

VO 700

60

k 0.005

15
g 9.8

m 3

Fig. Nº3 Lanzamiento de proyectil en 2D.

Problema Nº3

a. Determinar el modelo dinámico del sistema.

Parámetros:

• Velocidad inicial del lanzamiento: VO = 700

• Ángulo inicial con respecto a la horizontal: = 60

• Coeficiente de fricción: k = 0,005

• Aceleración de gravedad: g = 9,8

• Masa del proyectil: m = 3

Tenemos que considerar el efecto de la resistencia del aire en un lanzamiento parabólico, donde:

es la Fuerza del rozamiento del aire.

Luego por la segunda Ley de Newton se tiene:

−En el eje x:

(a1)

−En el eje y:

16
(a2)

Luego tenemos:

Que es representado por el siguiente sistema de ecuaciones:

b. Obtener de manera analítica las coordenadas X e Y

La ecuación (a1) es una ecuación diferencial, donde se obtiene la solución de la forma

. Así con las condiciones iniciales
y

obtenemos x.

La ecuación (a2) también es una ecuación diferencial, donde se obtiene la solución de la forma

. Así con las condiciones iniciales

y

obtenemos y.

En resumen:

17
c. Obtener una expresión aproximada para el tiempo de vuelo y la distancia recorrida.

De las ecuaciones del movimiento podemos obtener el tiempo de vuelo:

Despejando t de y(t)=0:

Aproximando:

Luego, desarrollando:

Lo que se aproximo fue:

Despreciando los términos

se llega a:

(a3)

18
Cuando
:

(a4)

Reemplazando (a4) en (a3):

Tiempo de vuelo aproximado.

Reemplazando los parámetros se llega a:

−Ahora podemos obtener la distancia recorrida

Ahora reemplazamos el tiempo de vuelo y los parámetros en esta ecuación, así tenemos:

d. Simule el sistema, obtenga gráficas para el desplazamiento [km] y velocidad [m/s], para el modelo
obtenido en a. y b. , con los parámetros dados en la tabla, y compare.

Se ocuparon los siguientes código para los graficos:

E3.m:

%Lunes 05/Octubre/2007

%Tarea 1 Modelación de Procesos: PROBLEMA N°3: LANZAMIENTO PROYECTIL

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

19
clear all;

close all;

clc;

global v0 theta k g m

%Parámetros del modelo:

v0=700;

theta=60*pi/180;

k=5e3;

%k=1e2

%k=2e2

%k=4e2

%k=8e2

g=9.8;

m=3;

%==========================================================================

%==========================================================================

%Condiciones Iniciales

x0=0;

vx0=v0*cos(theta);

y0=0;

vy0=v0*sin(theta);

%==========================================================================

%==========================================================================

%Gráfico:

[t,x]=ode45('P3_ec',[0 150],[x0 vx0 y0 vy0]);

figure(1)

20
plot(t,x(:,1)/1e3), axis on, grid on, title('x(t)')

xlabel('Tiempo'), ylabel('Posicion [Km]') %x

figure(2)

plot(t,x(:,2)), axis on, grid on, title('vx(t)')

xlabel('Tiempo'), ylabel('Velocidad [m/s]')%vx

figure(3)

plot(t,x(:,3)/1e3), axis on, grid on, title('y(t)')

xlabel('Tiempo'), ylabel('Posicion [Km]') %y

figure(4)

plot(t,x(:,4)), axis on, grid on, title('vy(t)')

xlabel('Tiempo'), ylabel('Velocidad [m/s]')%vy

y E3_ec.m:

function v_estado=E3_ec(t,x)

global v0 theta k g m

v0=700;

theta=60*pi/180;

k=0.005;

g=9.8;

m=3;

%==========================================================================

%==========================================================================

%Sistema de ecuaciones:

v_estado=zeros(4,1);

v_estado(1)=x(2);

v_estado(2)=−k*x(2);

v_estado(3)=x(4);

21
v_estado(4)=−k*x(4)−g;

Las graficas obtenidas son:

−Desplazamiento [Km]:

−Velocidad [m/s]:

e. Obtenga las trayectorias para los siguientes valores de k, con k= 0.005, 0.01, 0.02, 0.04 y 0.08.

Se ocupo el siguiente código:

E3e:

%Lunes 05/Octubre/2007

%Tarea 1 Modelación de Procesos: PROBLEMA N°3e: LANZAMIENTO PROYECTIL

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

clear all

close all

clc

%Parámetros del modelo:

v0 = 700;

theta=60*pi/180;

k = input('Introduzca coeficiente de fricción del aire: k= ');

g = 9.8;

m = 3;

v0y =v0*sin(theta);

v0x =v0*cos(theta);

t = 1:.2:120;

%==========================================================================

%Ecuaciones:

x = (v0x)*(1−exp(−k*t))/k;

y = −g*t/k+((g+k*v0y)/k^2)*(1−exp(−k*t));

22
%==========================================================================

%Gráfico:

figure(1),plot(x/1000,y/1000), axis on , grid on

xlabel('x [km]'), ylabel('y [km]') ,title('k=');

Y se obtuvieron las siguientes gráficas:

f. De la expresión aproximada en c. , para el tiempo de vuelo, cual es el rango de valor de k, para que la
aproximación sea válida.

Según las gráficas se puede apreciar que k tiene que ser del orden de 10−3 o menor para poder despreciar la
fuerza de roce del aire.

Si k> 10−3 no seria recomendable despreciar la fuerza de roce del aire en este caso.

Universidad de Concepción

Ingeniería Civil Electrónica

Universidad de Concepción

23