jose2997

Justin Drew Bieber primero que todo es este video de HALO que me costo mucho trabajo consegir ya que no sabia el nombre de la cancion link: http://www.youtube.com/watch?feature=endscreen&v=v4VNP4Q2NJw&NR=1 AL GUNA VEZ ME PREGUNTABA QUE SERIA LO MAS ENRETENIDO PARA Q VARIAS PERSONAS VIERAN VARIAS VECES EN INERNET YO EMPEZE A ODIAR ESTE JOVEN SIN SABER EN REALIDAD QUIEN ERA GRACIAS A LA CANTIDAD DE GIFS Y PUBLICIDAD Q DECIAN Q ESTE CHICO ERA GAY , Y PUBLICIDAD EN CONRA DE BIEBER PERO EMPEZE A OIR SU MUSICA Y ES BUENA SOBRE TODO SU VIDEO BABY PERO ESTE VIDEO NO ES EL UNICO EN LA CUENTA TOP 10 YA QUE SU VIDEO Never Say Never ft. Jaden Smith ACUPA EL PUESTO NUMERO DIEZ CON 290.917.758 vistas PERO LO Q ME MOTIVO A HACER ESTE POST SON LAS 700 000 000 VISITAS DE SU VIDEO BABY YA QUE NISIQUIERA LADY GAGA Y JLO LE LLEGAN A LOS TALONES CON LAS 470 000 000 VISITAS SE PUEDE DECIR QUE BIEBER LES LLEVA EL DOBLE NO ENCUENTRO TODAVIA PRUEBA VERDADES DE Q EL CHICO SEA GAY ASI Q NO DEVERIAN HACER MAS MOFA DE ESTO Nombre real Justin Drew Bieber Nacimiento 1 de marzo de 1994 (17 años) Origen London (Ontario), Bandera de Canadá Canadá Ocupación Cantante Información artística Género(s) Pop, R&B, teen pop Instrumento(s) Voz, batería, piano, guitarra, trompeta2 3 4 Período de actividad 2009-presente Discográfica(s) Island Records Artistas relacionados Usher Web Sitio web www.justinbiebermusic.com Señoras y señores 700 000 000 visitas en youtube Ficha Justin Bieber en IMDb link: http://www.youtube.com/watch?v=kffacxfA7G4&ob=av2e Señoras y señores 700 000 000 visitas en youtube sin mas palabra me despido para mi un hecho historico ya que LADY GAGA y JLO no alcanzan las 500 visitas seguire agregando cosas al post como un calculo matematico de cuando llegara a 800 000 000 visitas asi que siganlo si eres de los pocos fans de Justin Dejame puntos GRACias por ultimo la ñapa some doby to love en español alguna que amar.... link: http://www.youtube.com/watch?v=SOI4OF7iIr4

%----------------------------------------------------------------------------- % MIS CODIGOS PARA PRACTICAR IF FOR ELSE WHILE %----------------------------------------------------------------------------- (Al final puede verse el formato que se dio a esa práctica) %------------------------------------------------ function [L,U,p]=lupa(A,piv) %------------------------------------------------ % function [L,U,p] = lupa(A,piv) % Halla la factorización LU de A (que puede ser rectangular), % salvo que encuentre un pivote nulo. % Lo hace SIN pivotar si el argumento 'piv' es 0 (o no se da) % y PIVOTANDO si se da 'piv' distinto de 0 . En ese caso, % halla L,U tales que L*U = P*A , pero devuelve simplemente % un vector 'p' que da el orden de las filas en P*A. % Si size(A)=[m,n] , entonces size(L)=[m,m] , U es como A . if nargin<2, piv=0; end % si se ha llamado 'lupa(A)' [m,n] = size(A); q = 1:m; % filas del bloque restante p = q; % orden inicial de las filas %-------------------------------- loop de ELIMINACION: for k = 1:min(m-1,n) if piv % escoge la fila ip del PIVOTE: [a,ip] = max(abs(A(q,k))); ip = q(ip); p([k,ip]) = p([ip,k]); A([k,ip], = A([ip,k],; end if A(k,k)==0, error('pivote nulo'), end % !! q = k+1:m; col = k+1:n; % el bloque restante es A(q,col) A(q,k) = A(q,k)/A(k,k); % multiplicadores A(q,col) = A(q,col) - A(q,k)*A(k,col); end %-------------------------------- U = triu(A); L = A-U; L(m,m) = 0; % hay que agrandar L, si n < m L = L(:,1:m)+eye(m); %------------------------------------------------ function X = luaxb(A,B,piv) %------------------------------------------------ % function X = luaxb(A,B,piv) % Usa la factorización LU de A (llamando a 'lupa(A,piv)'), % para resolver el S.E.L. A*X = B . % B puede tener varias columnas, pero A debe ser cuadrada. % Lo hace SIN pivotar si el argumento 'piv' es 0 (o no se da) % y PIVOTANDO si se da 'piv' distinto de 0 . [m,n]=size(A); if m ~= n, error('A no es cuadrada'), elseif nargin<3, piv=0; end % ahora llama a la factorización LU [L,U,p]=lupa(A,piv); if prod(diag(U))==0, error('U no es regular'), end X=B(p,; % dará error si B tiene <m filas; % si tiene >m , las ignora for k=1:m-1 % ETAPA HACIA DELANTE ... q=k+1:m; q, = q, - L(q,k)*k,; end for k=m1 % ... Y ETAPA HACIA ATRAS k, = k,/U(k,k); q=1:k-1; q, = q, -U(q,k)*k,; end %------------------------------------------------ function [Q,R] = qra(A,met) %------------------------------------------------ % function [Q,R]=qra(A,met) % % Produce los factores QR de A usando Gram-Schmidt clásico si met=0, % y en caso contrario, Gram-Schmidt modificado. % Admite A de cualquier forma y rango, y devuelve factores "reducidos": % con r columnas en Q y filas en R, r=rango(A). if nargin<2, met=1; end [m,n] = size(A); Q=zeros(m,0); % para empezar en k=1 sin que " Q'*v " dé error R=0*A; K=[]; % K = Ãndices de filas R y cols Q que se conservarán if met==0 %------------------------------------ Gram-Schmidt clásico for k=1:n v = A(:,k); j=1:k-1; R(j,k) = Q'*v; v = v -Q*R(j,k); d = norm(v); if d > 0 Q(:,k) = v/d; R(k,k) = d; K=[K,k]; end end else %------------------------------------ Gram-Schmidt modificado for k=1:n q = A(:,k); d = norm(q); if d > 0 q = q/d; Q(:,k) = q; j= (k+1):n; p = q'*A(:,j); R(k,k:n) = [d,p]; A(:,j) = A(:,j) -q*p; K=[K,k]; end end %------------------------------------ end % retira lineas inútiles: Q = Q(:,K); R = R(K,; %------------------------------------------------ function [W,R,Q]=house(A) %------------------------------------------------ % % function [W,R,Q]=house(A) % Halla la factorización QR de A mediante reflexiones de Householder, % y devuelve el factor R y los vectores unitarios W_j que realizan % las reflexiones H_j , como columnas de la matriz W . % Si se la llama con el tercer argumento de salida, devuelve en él % el factor ortogonal Q , producido a partir de la W . [m,n] = size(A); %------------------------------------ Householder K=[]; % K = índices de las W_k no nulas for k= 1:min(m-1,n) % es el número de columnas a "limpiar" i=k:m; % sólo hay que operar con estas filas v = A(i,k); d = norm(v); if d > 0 % si fuese v=0, nada que hacer esta vez K=[K,k]; sign = 1-2*(v(1)>0); % signo = -1 si v(1)>0 A(k,k) = sign*d; v(1) = v(1) - A(k,k); v = v/norm(v); W(i,k) = v; j= (k+1):n; % columnas en las que opera H_k : A(i,j) = A(i,j) - 2*v*(v'*A(i,j)); end end %------------------------------------ W=W(:,K); % limpia sobrantes en W y R R=triu(A); if nargout>2 % si se ha pedido, calcula Q ... Q = eye(m); % ... aplicando a I las H_k ... for j = size(W,2)1 % ... en orden inverso v = W(:,j); Q = Q - 2*v*(v'*Q); end end %------------------------------------------------ function X = qraxb(A,B,met) %------------------------------------------------ % function X = qraxb(A,B,met) % % Usa la factorización QR de A para resolver el S.E.L. A*X = B , % llamando % a 'qra(A,met)' si se da met = 0 ó 1 % a 'house(A)' si no se da ese argumento de entrada. % % Si A es rectangular, pero con rango = no.de cols , calcula por % ese método la solución LS . if nargin<3 %------------------------------------ [W,R]=house(A); for j = 1:size(W,2) % v = W(:,j); B = B - 2*v*(v'*B); end X = RB; elseif met==0 %------------------------------------ Gram-Schmidt clásico [Q,R]=qra(A,met); B = Q'*B; X = RB; else %------------------------------------ Gram-Schmidt modificado [Q,R]=qra(A); B = Q'*B; X = RB; %------------------------------------ end %------------------------------------------------ function [X,res,normrk] = iteraxb(A,B,w,tol) %------------------------------------------------ % function [X,res] = iteraxb(A,B,w) % % Usa un método iterativo para el S.E.L. A*X = B : % Jacobi, si se da w = 0, % SOR, si se da un valor 0 < w < 2 ; en particular, % Gauss-Seidel si se da w = 1 . % Devuelve la aproximación X y su correspondiente res = B - A*X . % % Condiciones de parada: % que el 'res' crezca a un tamaño muy superior al de B ; % que descienda hasta una fracción 'tol' de ese tamaño, que % puede darse como otro argumento de entrada (por defecto, 1/10^6). % Si se la llama con un 3er argumento de salida, devuelve el % historial de la convergencia: las sucesivas norm(r_k)/norm(B). if nargin < 4, tol=10^-6; end m = size(B,1); if any( size(A) ~= m ), error('datos inadecuados'), end tol = tol* norm(B); % tamaño de parada para el residual it = 1; M = 100; % contador y máximo de iteraciones if w==0 S = diag(diag(A)); %----------------------- Jacobi else S = tril(A); %------------------------ SOR end X = 0*B; % valor inicial x = 0 res = B; % residual inicial normrk = norm(res); while normrk(it) > tol & it < M if w==0, X = X + (Sres); else, X = X + w* (Sres); end res = B - A*X; it = it+1; normrk(it) = norm(res); end if nargout > 2, normrk = normrk/norm(B); end %------------------------------------------------ SESION 3 DE LAB (J 25- L 29/10) %-------------------------------------------------------------------- % EXPERIMENTOS CON 'lu' % Producir una A a partir de sus factores LU, % de tal modo que NO HAYA QUE PIVOTAR ?? m=5; % m = tamaño escogido format bank d=5; % Cortaremos'LU' para producir L,U , LU = ceil(d*rand(m))/d % ... y no queremos un |l_ik|>1 , % ... ni ceros en la diagonal de U. A = (tril(LU,-1)+eye(m))*(triu(LU)) [L,U] = lu(A) % y Matlab recupera los factores L,U format short %-------------------------------------------------------------------- % RESOLVIENDO UN S.E.L. Ax=b con 'Ab'... % Producimos b=A_1+A_2 : z=zeros(m,1); z(1:2)=1; b = A*z; x = Ab % y Matlab recupera x ... err = x-z % ... con errores ... res = b-A*x % consistentes con los k(L), k(U), k(A): cnLUA = [cond(L),cond(U),cond(A)] % (!!) Si se toma d=4 al producir A (bloque anterior) los errores ... ?? %-------------------------------------------------------------------- % ... Y CON LOS METODOS ITERATIVOS CLASICOS m=5; A=rand(m)/m+eye(m) % Ahora necesitamos A para las que cn = cond(A) % ... el método converja z=zeros(m,1); z(1:2)=1; b = A*z; % Jacobi: S = A.*eye(m); T = S-A; normaDG = norm(ST) x = 0*b; for it=1:10 xj = S(T*x+b); err(it) = norm(xj-x), res(it) = norm(b-A*xj) x=xj, pause end % Gauss-Seidel: S = triu(A); T = S-A; normaDG = norm(ST) x = 0*b; for it=1:10 xj = S(T*x+b); errG(it) = norm(xj-x); resG(it) = norm(b-A*xj); x=xj; end error = norm(x-z), residual = norm(b-A*x) it=1:10; semilogy(it,err,'*k',it,res,'.k',it,errG,'*r',it,resG,'.r') % Comparemos el rendimiento de ambos, m=400; % ... con m "grande" A=rand(m)/m+eye(m); cn = cond(A) z=zeros(m,1); z(1:2)=1; b = A*z; % El método DIRECTO: (Matlab usa LU) tic, x=Ab; toc % ... frente a Gauss-Seidel: S = triu(A); T = S-A; normaDG = norm(ST) tic for it=1:10, x = S(T*x+b); end toc error = norm(x-z), residual = norm(b-A*x) %-------------------------------------------------------------------- % PERDIDA INESPERADA DE RANGO: las L de LU son MUY especiales!! for m=40:20:120 for it = 1:20 i = it+m-40; A = rand(m); % (!!) probar con A = randn(m); [L,U,P] = lu(A); kn(i, = [cond(L),cond(U),cond(A)]; ra = [rank(L),rank(U),rank(A)]; if min(ra)<m, ra, end LA = tril(A); K(i) = cond(LA); r(i) = rank(LA); end end i = 1:i; m = 40+fix((i-1)/20)*20; % figure subplot(2,1,1), semilogy(i,kn,i,m,'k') ylabel('rango y cond de las L, U, A') subplot(2,1,2), semilogy(r,K,'r.'), grid xlabel('rango de las "L" recortadas: ligeramente deficitario') ylabel('sus "cond"') %------------------------------------------------ SESION 5 DE LAB (L 12, J 15/11). %-------------------------------------------------------------------- % A) fase de prueba del código con mis propias versiones de las 6 funciones: % pruebas con matrices A pequeñas (que recomiendo repetir con las de cada uno), % para ver simplemente si funcionan como es debido. %-------------------------------------------------------------------- m=4; A=rand(m); clc %----------------------- Factorización LU : [L0,U0,P] = lu(A) % con la de Matlab [L,U,p] = lupa(A,1) % con la nuestra bkerr = A(p, - L*U help lupa, pause, disp('... y sin pivotar:') [L,U] = lupa(A); bkerr = A - L*U , pause, clc %----------------------- Factorización QR : [Q0,R0] = qr(A) % la de Matlab [Q,R] = qra(A) % la que probamos: Gr-Sch modificado, ... disp(' OJO A LOS SIGNOS !!'), bkerr = A - Q*R pause, disp('... ahora, Gr-Sch clásico:') [Q,R] = qra(A,0) bkerr = A - Q*R , help qra pause, disp('... y ahora, Householder:') [W,R,Q] = house(A) bkerr = A - Q*R , help house, pause, clc %----------------------- Resolver Ax=b : x0=zeros(m,1); x0(1:2)=1, b=A*x0; xM = Ab; % con Matlab errMatlab = xM - x0 x = luaxb(A,b,1); % con LU errpiv = x - x0, help luaxb, pause x = luaxb(A,b,0); errNOpiv = x - x0, pause x = qraxb(A,b,1); % llamando a 'qra' errGSch = x - x0, help qraxb, pause x = qraxb(A,b); % llamando a 'house' errHouse = x - x0, pause, clc %----------------------- M»todos iterativos A=A+sqrt(m)*eye(m); b=A*x0; format short e, clc [X,res,hist] = iteraxb(A,b,0), help iteraxb, pause tol = 10^-10 w=1, [X,res,hist] = iteraxb(A,b,w,tol); hist S = tril(A); DG = S(S-A); ep=10^-6; plot(eig(DG)+ep*i,'*r'), axis equal, hold on format short, rho = max(abs(eig(DG))) t = 0:0.01:1; z = exp(j*2*pi*t); plot([z,rho*z]), hold off %-------------------------------------------------------------------- % B) fase de exploración (comienzo): %-------------------------------------------------------------------- %----------------------- Métodos iterativos m=40; d=5; U=triu(ones(m),1); U=U^d; U=U+U'; P=1-(U>0); A=randn(m)+sqrt(m)*eye(m); A=A.*P; spy(A), pause x0=zeros(m,1); x0(1:2)=1; b=A*x0; tol = 10^-10 w=1, [X,res,hist] = iteraxb(A,b,w,tol); hist, % rhow = max(abs(1+w*(eig(DG)-1))) S = tril(A); DG = S(S-A); format short, rho = max(abs(eig(DG))), ep=10^-6; plot(eig(DG)+ep*i,'*r'), axis equal, hold on t = 0:0.01:1; z = exp(j*2*pi*t); plot([z,rho*z]), hold off %----------------------- Factorización QR : m=100; A=rand(m); [Q0,R0]=qr(A); % R0=R0*diag(sign(diag(R0))); A=Q0*R0; % R0=triu(A)*diag(sign(diag(A))); A=Q0*R0; kA = cond(A) disp('ERRORES:') disp('norm(Q*R-A)/norm(A) , cond(Q)-1 , norm(Q-Q0) , norm(R-R0)/norm(R0)') format short e disp(' -------------------- Gr-Sch CLASICO:') [Q,R] = qra(A,0); errores0 = [norm(Q*R-A)/norm(A), cond(Q)-1, norm(Q-Q0), norm(R-R0)/norm(R0)]; disp(errores0), pause disp(' -------------------- modificado:') [Q,R] = qra(A); errores1 = [norm(Q*R-A)/norm(A), cond(Q)-1, norm(Q-Q0), norm(R-R0)/norm(R0)]; disp(errores1), pause disp(' -------------------- Householder:') [W,R,Q] = house(A); erroresH = [norm(Q*R-A)/norm(A), cond(Q)-1, norm(Q-Q0), norm(R-R0)/norm(R0)]; disp(erroresH), pause disp(' -------------------- qr de Matlab:') [Q,R] = qr(A); erroresM = [norm(Q*R-A)/norm(A), cond(Q)-1, norm(Q-Q0), norm(R-R0)/norm(R0)]; disp(erroresM), format short %----------------------- Resolver Ax=b : x0=zeros(m,1); x0(1:2)=1; B=A*x0; disp(' -------------------- Householder:') X = qraxb(A,B); errH = norm(X-x0), pause disp(' -------------------- Gr-Sch CLASICO:') X = qraxb(A,B,0); err0 = norm(X-x0), pause disp(' -------------------- Modificado:') X = qraxb(A,B,1); err1 = norm(X-x0), pause disp(' -------------------- LU:') X = luaxb(A,B); err0 = norm(X-x0), pause disp(' -------------------- LU con pivotaje:') X = luaxb(A,B,1); err1 = norm(X-x0), pause disp(' -------------------- Ab de Matlab:') errMatlab = norm((AB)-x0) %-------------------------------------------------------------------- % PLAZO Y FORMA PARA LA PRACTICA 2: %-------------------------------------------------------------------- Deberá estar entregada al mediodía (12:00) del V 23/11, dejando cada equipo EN LA CARPETA 'practica2**', QUE SE LES VA A CREAR PARA ELLO, los documentos siguientes: * 6 M-files con las funciones cuyos nombres y uso se especifica abajo; * un documento de texto que explique (con toda la concisión que lo complejo del asunto permite) las pruebas que han hecho con el código, comentando los resultados que consideren más relevantes o hayan llamado su atención, en relación con los siguientes aspectos: - costes: tiempos de ejecución, comparados entre sí y/o con los de rutinas equivalentes de Matlab ('lu(A), qr(A), Ab'); - precisión y estabilidad: errores "back-/forward" y su relación con los 'cond' de las diversas matrices; - y la evolución de ambas cosas al crecer el tamaño de las matrices. En la sesión de esta semana próxima ofreceré ejemplos de ese tipo de indagación (usando mis M-files); naturalmente, se pueden escribir varios volúmenes sobre algo así (más aún con todos los métodos que abarcamos en la Práctica), de modo que una selección prudente será parte importante del mérito del trabajo. Las 6 funciones deben tener los nombres, forma y comportamiento que se explica abajo. Las especificaciones precedidas por [*] son "OPTATIVAS": se debe conseguir como mínimo que los códigos funcionen para una matriz cuadrada y regular, y añadir esas otras opciones sólo si sobra tiempo y ganas. Para verificar que funcionan como es debido, se recomienda usar datos de tamaño pequeño (n=4, por ejemplo) producidos con 'rand' y/o 'randn' para poder visualizar los resultados, además de compararlos con los que calcula Matlab. Esa 'fase A' de las pruebas NO DEBE formar parte del "informe" que acompañe a los códigos, en el que debe comentarse sólo lo que pasa al probar con datos más "exigentes". ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! NO OLVIDAR incluir el NUMERO DE EQUIPO Y NOMBRES DE SUS COMPONENTES en las lineas iniciales de comentario de CADA FUNCION, que deben seguir inmediatamente a la linea de declaración de la misma, para que Matlab las reconozaca y use como su 'help'. ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! Y naturalmente, redactar esos comentarios pensando en que sirvan para eso: copiar los que pongo a continuación y añadir, además de los DATOS citados, lo que proceda explicar si se han programado más opciones de entrada o salida. %-------------------------------------------------------------------- function [L,U,p]=lupa(A,piv) % % Halla la factorización LU de A (cuadrada), % pero para y devuelve 'error' si encuentra un pivote nulo. % Lo hace SIN pivotar si el argumento 'piv' es 0 (o no se da) % y PIVOTANDO si se da 'piv' distinto de 0 . En ese caso, % halla L,U tales que L*U = P*A , pero devuelve simplemente % un vector 'p' que da el orden de las filas en P*A. %-------------------------------------------------------------------- function X = luaxb(A,B,piv) % % Usa la factorización LU de A (llamando a 'lupa(A,piv)'), % para resolver el S.E.L. A*X = B . % Por lo tanto, lo hace SIN pivotar si el argumento 'piv' es 0 % (o no se da) y PIVOTANDO si se da 'piv' distinto de 0 , % y devuelve 'error' si 'lupa' encuentra algún pivote nulo. % % B puede tener varias columnas, pero A debe ser cuadrada. %-------------------------------------------------------------------- function [Q,R]=qra(A,met) % % Produce los factores QR de A usando Gram-Schmidt clásico si met=0, % y en caso contrario, Gram-Schmidt modificado. % [*] Admite A de cualquier forma y rango, y devuelve factores "reducidos": % con r columnas en Q y filas en R, r=rango(A). %-------------------------------------------------------------------- function [W,R]=house(A) % % Halla la factorización QR de A mediante reflexiones de Householder, % y devuelve el factor R y los vectores unitarios W_j que realizan % las reflexiones, como columnas de la matriz W . % [*] Si se la llama con un tercer argumento de salida, devuelve en él % el factor ortogonal Q , producido a partir de la W . %-------------------------------------------------------------------- function X = qraxb(A,B,met) % % Usa la factorización QR de A para resolver el S.E.L. A*X = B , % llamando % a 'qra(A,met)' si se da met = 0 ó 1 % a 'house(A)' si no se da ese argumento de entrada. % [*] Si A es rectangular, pero con rango = no.de cols , calcula por % ese método la solución LS . %-------------------------------------------------------------------- function [X,res] = iteraxb(A,B,w) % % Usa un método iterativo para el S.E.L. A*X = B : % Jacobi, si se da w = 0, % SOR, si se da un valor 0 < w < 2 ; en particular, % Gauss-Seidel si se da w = 1 . % Devuelve la aproximación X y su correspondiente res = B - A*X . % % Condiciones de parada: % que el 'res' crezca a un tamaño muy superior al de B ; % que descienda hasta una fracción fijada de ese tamaño. % [*] Esa fracción puede darse como un 4o. argumento de entrada 'tol'. % Si se la llama con un 3er argumento de salida, devuelve el historial % de la convergencia: la norma de los sucesivos residuales /norm(B) .

% CONTROL MANUAL Y AUTOMÁTICO DE PÉNDULO INVERTIDO % % Este script ejecuta la simulación de un péndulo que puede ser controlado manualmente (accionando % una barra deslizante o "slider" o de forma automática, mediante lazos de control PID. El control % automático incorpora un doble lazo incluyendo un lazo con control del ángulo del péndulo (para que se sitúe en % posición vertical) y otro lazo más lento para hacer que la posición de la base se sitúe en un punto dado % en el estado de equilibrio. % % El script puede ser utilizado tanto como un "equipo de prácticas", en el que pueden % estudiarse diversos conceptos de control, así como servir de base para que el profesor (o el % alumno) elaboren otros ejemplos interactivos totalmente distintos mediante el procedimiento % de "cortar y pegar". % % % Fecha: 2006-11-02 % Autor: Ignacio Díaz % Area de Ingeniería de Sistemas y Automática % Universidad de Oviedo clear; close all; clc; disp('Instrucciones:'); disp('- Pulsar ''0'' para desactivar el sistema de control'); disp('- Pulsar otra tecla para restaurar el sistema de control'); disp(' '); disp('Control manual: actuar sobre la barra deslizante para modificar la base del péndulo.'); pause(2); % PARÁMETROS DEL PÉNDULO l = 2; % Longitud del péndulo m = 1; % Masa del péndulo J = m*l^2; % Momento de inercia referido al eje B = 1; % Coeficiente de fricción g = 10; % Aceleración de la gravedad % ESTADO INICIAL DEL PÉNDULO x = [pi-0.1;0]; % Para que se vea el efecto del control, empezamos % con el péndulo casi vertical (theta = pi +/- "algo" % DEFINICIÓN DE UN "SLIDER" PARA CONTROLAR MANUALMENTE EL PÉNDULO f = figure(1); set(f,'pos',[200,200,700 700]); h = uicontrol('style','slider','pos',[20 20 680 20],'min',-4,'max',4); Tm = 0.01; % Período de muestreo e0 = zeros(2,1); % Condiciones iniciales del control de ángulo ep0 = zeros(2,1); % Condiciones iniciales del control de posición x0 = [0;0]; % Condiciones iniciales del péndulo a0 = [0;0]; xmin = -2; xmax = +2; y = x(1); % PARÁMETROS DEL CONTROL PID DEL ÁNGULO Kp = 10; Ki = 100; Kd = .1; s = tf('s'); Cth = c2d(Kp + Ki/s + Kd*s/(0.01*s+1),Tm,'tustin'); % PARÁMETROS DEL CONTROL PID DE LA POSICIÓN DE LA BASE Cpos = c2d(.04 + 0.0001/s + s*0.0001/(0.01*s+1),Tm,'tustin'); pos = 0; % Valor inicial de la posición de la base del péndulo k = 2; % Empezamos en k=2 para tener acceso al menos a dos muestras anteriores while 1, k = k + 1; % BUFFER CON LA POSICIÓN (3) TIENE LA POSICIÓN ACTUAL) % k) = get(h,'Value'); e(k) = pi-y; % RESTRINGIMOS EL VALOR DEL ERROR AL INTERVALO [-pi,pi] e(k) = mod(e(k)+pi,2*pi)-pi; % CONTROL DE LA POSICIÓN DE LA BASE DEL PÉNDULO ep = 0-pos; [dpos,ep0] = filter(Cpos.num{1},Cpos.den{1},ep,ep0); e(k) = e(k) - dpos; % SELECCIÓN DE CONTROL MANUAL / CONTROL AUTOMÁTICO tecla = get(f,'currentchar'); switch tecla case '0' pos = get(h,'value'); otherwise [pos,e0] = filter(Cth.num{1},Cth.den{1},e(k),e0); end % Suavizado de la aceleración (muy conveniente, porque el movimiento del % objeto "slider" con un ratón se produce a saltos, dando lugar a segundas % derivadas muy elevadas) [k),x0] = filter(.01,poly([.9 .9]),pos,x0); [a,a0] = filter((1/Tm^2)*[1 -2 1],[1 0 0],k),a0); A(k) = a; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % ECUACIONES EN ESPACIO DE ESTADOS (NO LINEALES) DEL PÉNDULO u = -a; % Asignamos la entrada % Ecuación de estados x(1) = x(1) + Tm*x(2); x(2) = x(2) + Tm*(1/J*(-B*x(2)-m*g*l*sin(x(1))+m*u*l*cos(x(1)))); % Ecuación de salida y = x(1); th = y; % Asignamos la salida %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % REPRESENTACIÓN GRÁFICA DE LA SIMULACIÓN figure(1); plot(k),0,'.'); hold on; p1 = k); p2 = k)+l*exp(j*(th-pi/2)); line(real([p1,p2]),imag([p1,p2])); plot(real(p2),imag(p2),'.','markersize',40); hold off; % Sugerencia: pueden dibujarse también otras flechas indicando en tiempo real las fuerzas reales % o de inercia que actúan en cada elemento del sistema % Centrado automático de la perspectiva sobre el objeto de control if k)>xmax-1 xmin = xmin + 0.1; xmax = xmax + 0.1; elseif k)<xmin+1 xmin = xmin - 0.1; xmax = xmax - 0.1; end grid on; axis([xmin-3 xmax+3 -5 5]); % Refresco de la imagen drawnow; end