Minicurso
MATLAB
Rodrigo Cardim
Laboratório de Pesquisa em Controle
DEE-FEIS-UNESP
24-25/05/2010
1
1 - Conceitos Básicos
Material auxiliar: Introdução ao MATLAB
Prof. Edvaldo Assunção
Prof. Ricardo Tokio Higuti
DEE-FEIS-UNESP
2
Introdução
MATrix LABoratory
Inicialmente escrito em FORTRAN
Novo MATLAB escrito em C
Várias plataformas
Win, Unix, Linux, Macintosh
3
Versões
4
Características:
Álgebra matricial
Versatilidade: O usuário cria novas
ferramentas
Programação com macro funções
ToolBoxes específicos
Vários livros baseados em Matlab
Atualmente na versão 7.10 (R2010a)
5
Janela Principal do Matlab
6
Comandos
= Comando de atribuição
a=10
[ ] Delimita elementos de matrizes e vetores
x=[0 1 2 3] y=[0 1; 2 3]
% Comentário
%comentário
help Tópicos de ajuda
help cos
: Incremento
7
Exemplo:
Digite as seguintes linhas de comando
>> % definição de um vetor
>> b1=[1 2 3 4 5 6 7 8 9]
>> b2=[1;2;3;4;5;6;7;8;9]
>> %definição de uma matriz
>> c=[1 2 3;4 5 6;7 8 9]
>> c=[c;[10 11 12]]
>> c(2,2)=0
>> % Criação de vetores com incremento
>> x=1:2:9
>> x=0:pi/3:pi;
>> y=sin(x)
>> c(:,1)
>> c(1:2,1)
>> clear clc .m 8
Comando format
format short ponto fixo com 4 casas decimais (o padrão) 35.8333 (default)
format long ponto fixo com 14 casas decimais 35.83333333333334
format short e notação científica com 4 casas decimais 3.5833e+01
format long e notação científica com 15 casas decimais 35.83333333333334e+01
format hex hexadecimal 4041eaaaaaaaaaab
format bank duas casa decimais 35.83
format rat aproximação racional 215/6
format + positivo, negativo ou zero +
realmin 2.2251e-308
realmax 1.7977e+308
9
Operações com Matrizes e Vetores
+ Adição
- Subtração
* Multiplicação
/ Divisão à direita
\ Divisão à esquerda
^ Exponenciação
’ Transposta
inv( ) Inversão matricial
10
Transposta, adição e subtração
>> a=[1 2 3;4 5 6;7 8 9];
>> b=a’
>> c=a+b
>> c=a-b
Multiplicação e adição com escalar
>> x=[-1 0 2];
>> y=[-2 -1 1]’;
>> x*y
>> c=x+2
11
Inversão e Divisão
>> a=[1 0 2;0 3 4;5 6 0];
>> b=inv(a)*a
>> c=b/a % c=b*inv(a)
>> c=b\a % c=inv(b)*a
12
Resolvendo Sistemas Lineares
x1 + 2x2 + 0x3 = 5 1 2 0 x1 5
1 5 3 . x 0
-x1 +5x2 - 3x3 = 0 2
4x1 - 2x2 + x3 = 3 4 2 1 x3 3
-1
A . X =B Solução : X = A . B
Resolvendo com o Matlab
>> A=[1 2 0;-1 5 -3;4 -2 1];
>> B=[5; 0; 3];
>> X=A\B
ou
>> X=inv(A)*B 13
Operação Elemento a Elemento – Matriz e Vetor
.* Multiplicação
./ Divisão à direita
.\ Divisão à esquerda
.^ Exponenciação
>> x=[1 -2 3];
>> y=[4 3 2];
>> z=x.*y
>> z=x.^y
>> y.^2
14
Operação com Números Complexos
>> z=3+4*i % i=sqrt(-1)
>> a=[1 2;3 4]+i*[5 6;7 8]
>> Mz=abs(z)
>> Az=angle(z)
15
Utilitários para Matrizes
>> a=eye(3) % Matriz identidade
>> a=zeros(4) % Matriz nula
>> a=ones(3) % Matriz unitária
>> a=rand(2,3) % Matriz com n. aleatórios
>> a=[2 0 0;0 3 0;0 0 -1];
>> d=det(a) % Determinante da matriz a
16
Autovalores e Autovetores
[Link] = li .xi
>> a=[1 0 0;0 -2 0;0 0 3];
>> l=eig(a)
>> a=[1 0 2;0 -2 0;0 1 3];
>> [x,l]=eig(a)
17
Traçando Gráficos
Gráfico do tipo y(t) x t
>> t=0:0.07:6*pi;
>> y=sin(t);
>> plot(t,y,’k’)
>> xlabel(‘tempo [s]’)
>> ylabel(‘sen(t)’)
.m
18
Duas ou mais curvas do tipo y(t) x t
>> z=cos(t);
>> plot(t,y,’b’,t,z,’r-.’)
>> title(‘Funções trigonométricas’)
>> xlabel(‘Tempo [s]’)
>> ylabel(‘Sen(t) e Cos(t)’)
>> text(3,0.6,’Seno’)
>> text(2.2,-0.5,’Cosseno’)
.m
19
>> t=0:pi/20:2*pi;
>> plot(t,sin(t),'-.r*')
>> hold on
>> plot(t,sin(t-pi/2),'--mo*')
>> plot(t,sin(t-pi/2),'--mo')
>> plot(t,sin(t-pi),':bs')
>> hold off
20
Gráfico Tridimensional
>> x=-8:0.5:8;
>> y=x’;
>> X=ones(size(y))*x;
>> Y=y*ones(size(x));
>> R=sqrt(X.^2+Y.^2)+eps;
>> Z=sin(R)./R;
>> surf(X,Y,Z);
>> xlabel(‘eixo X’);
>> ylabel(‘eixo Y’);
>> zlabel(‘eixo Z’)
>> title(‘Chapéu Mexicano’);
>> grid on;
21
>> x=-pi:0.2:pi;
>> y=-pi:0.2:pi;
>> [X,Y]=meshgrid(x,y);
>> Z=X.*sin(2*Y);
>> surf(X,Y,Z)
.m
22
Comando comet
.m
23
Comando comet3
.m
24
Aquisição de Dados com Osciloscópios
Tektronix
Obtendo dados do canal 1
>> [t,v]=curva(1);
>> tmin=min(t);
>> figure(1); plot(t-tmin,v)
>> xlabel(‘t (s)’)
>> ylabel(‘volts’)
>> title(‘Curva - Canal 1’)
>> grid on
25
Salvando os dados
>> save DadosCanal1 % gera um arquivo .mat
Obtendo dados do canal 1 e 2
>> [t,v]=curva([1,2]);
>> tmin=min(t);
>> figure(2); plot(t-tmin,v)
>> xlabel(‘t (s)’)
>> ylabel(‘volts’)
>> title(‘Curva - Canal 1 e 2’)
>> grid on
>> save DadosCanal12
26
Carregando dados salvos
>> load DadosCanal1 % Arquivo .mat
>> figure(1); plot(t-tmin,v)
>> xlabel(‘t (s)’)
>> ylabel(‘volts’)
>> title(‘Curva - Canal 1’)
>> grid on
.m
27
Comandos de Controle de Fluxo
O comando “for”
Formato:
for i=expressão
comandos;
end
28
Exemplo: arquivo teste1.m .m
n=3;m=3;
for i=1:m
for j=1:n
a(i,j)=i+j;
end
end
s=sprintf(‘\nMatriz A:a(i,j)=i+j\n’);
disp(s);disp(a)
>> teste1 % execução do programa
29
O comando “while”
Formato:
while condição
comandos;
end
30
Exemplo: arquivo teste2.m .m
n=1;
while n<=23
n=n+1;
end
disp(sprintf(‘\n n final: %d’,n));
No Matlab digite
>> teste2
31
O comando “if”
Formato:
if condição
comandos1;
else
comandos2;
end
32
Exemplo: arquivo teste3.m .m
%Este programa determina se o num. n é par ou ímpar
for n=1:4
resto=rem(n,2);
if resto==0
disp(sprintf(‘\n %d é par\n’,n));
else
disp(sprintf(‘\n %d é ímpar\n’,n));
end
end
33
Criando Subrotinas
Deve-se criar um arquivo com extensão .m
A subrotina deve ter o mesmo nome do
arquivo que a define
34
Arquivo com a função media.m .m
>> function x=media(u)
>> % Esta função calcula a média dos elementos de u
>> x=sum(u)/length(u);
Arquivo teste4.m .m
>> v=1:1:10;
>> m=media(v);
>> disp(sprintf(‘\n A média de 1 a 10 é: %4.2f’, m));
No MATLAB digite
>> teste4
35
Saindo do Matlab
>> quit
ou
>> exit
36
2 – Matlab e Simulink
(Exemplos)
37
Ambiente Simulink
38
Exemplo – Controle de um Motor DC
Matlab/Simulink
39
Simulação 1 – Sistema não forçado
.m .mdl
40
Simulação 2 – Controle Proporcional
.mdl
41
Simulação 3 – Controle PI
.mdl
42
43
3 – Simulações Utilizando a
Função ODE45 do MATLAB
44
Simulações Utilizando a Função ODE45 do
MATLAB
As EDOs são indispensáveis para o modelamento de muitos
sistemas físicos;
Descrevem como a taxa de variação das variáveis de um
sistema é influenciada pelas próprias variáveis do sistema e
por estímulos externos;
Naqueles casos em que as equações não podem ser
resolvidas prontamente de forma analítica, é conveniente
resolvê-las numericamente.
45
Exemplo
x (1 x 2 ) x x 0
Assim como em toda abordagem numérica para resolução de
equações diferenciais, as equações diferenciais de ordem
superior devem ser rescritas em termos de um conjunto
equivalente de equações diferenciais de primeira ordem.
y1 x y1 y2
y2 x y 2 (1 y12 ) y2 y1
Para resolver equações como essas, pode utilizar a função
ODE45 do Matlab!
46
Sintaxe
[T,Y] = ODE45(‘yprime’, Tinicial, Tfinal, Y0)
Com esta sintaxe integra-se um sistema de equações
diferenciais ordinárias descrito por um arquivo .M (yprime.m),
sobre um intervalo de tempo entre "Tinicial" e "Tfinal" e com
condições iniciais fornecidas em “Y0”.
47
Exemplo – Simulações com o Pêndulo Invertido
48
Modelos
sen( )(l cos( )m Mg mg ) u cos( )
l (m cos 2 ( ) M m)
Não-Linear
m sen( )( g cos( ) l ) u
x
m cos 2 ( ) M m
0 1 0 0 0
x1 x1 M m x1 1
x x g 0 0 0
Linear x 2 x
2 Ml . x2 Ml .u
x x3 x3 0 0 0 1 x3 0
x4 m g
0 0 0 4
x 1
x x4 M M
49
Controle
u Kx
Considerando:
- comprimento da haste: l=0,5m
- massa da bola: m=0,1Kg
- massa do carro: M=2,0Kg
- aceleração da gravidade: g=9,8m/s2
0 1 0 0 0
20,58 0 0 0 1
A e B
0 0 0 1 0
0, 49 0 0 0 0,5
50
Pólos desejados:
p1 2 2i, p2 2 2i, p3 5 2, p4 5 2
Utilize o comando acker do Matlab
0 1 0 0 0
20.58 0 0 0 1
A e B
0 0 0 1 0
0.49 0 0 0 0.5
>> K=acker(A,B,P)
K = -131.8861 -28.5152 -30.6122 -23.0892
Digite: eig(A-B*K)
51
Simulação – Modelo Não-Linear
%------------ Simulação do modelo Não-Linear a uma dada condição inicial -------------
%----------------------------------------- Arquivo 1 -----------------------------------------------
function dx=teste42(t,x)
%Esta função representa matematicamente o pêndulo invertido
%(sem linearização) envolvendo o sinal de controle u, sendo:
%x(1)=theta x(2)=theta' x(3)=x x(4)=x' u=-K1x(1)-K2x(2)-K3x(3)-K4x(4)
%considerando: .m
M=2; m=0.1; l=0.5; g=9.8
K1=-131.8861; K2=-28.5152; K3=-30.6122; K4=-23.0892;
dx=zeros(4,1)
dx(1)=x(2)
dx(2)=(sin(x(1))*(l*cos(x(1))*m*x(2)-g*M-g*m)+(cos(x(1))*(-K1*x(1)-K2*x(2) - K3*x(3)-K4*x(4))))/l/(-M-m+m*cos(x(1))^2)
dx(3)=x(4)
dx(4)=(m*sin(x(1))*(g*cos(x(1))-l*x(2))+K1*x(1)+K2*x(2)+K3*x(3)+K4*x(4))/ (-M-m+m*cos(x(1))^2)
%----------------------------------------- Arquivo 2 ------------------------------------------------
tempo=[0 4] % ou tempo=0:0.001:4 (fixa o passo);
cond_inicial=[0.6 0 0 0] %cond_inicial=[theta(rad) theta' x x']
[T,X]=ode45('teste42',tempo,cond_inicial)
subplot(2,2,1); plot(T,X(:,1))
xlabel('t(s)'); ylabel('theta(rad)'); title ('ângulo') .m
subplot(2,2,2); plot(T,X(:,2))
xlabel('t(s)'); ylabel('dtheta(rad/s)'); title ('Velocidade angular')
subplot(2,2,3); plot(T,X(:,3))
xlabel('t(s)'); ylabel('x(m)'); title ('posição')
subplot(2,2,4); plot(T,X(:,4))
xlabel('t(s)'); ylabel('dx(m/s)'); title ('Velocidade')
52
Simulação – Modelo Não-Linear
â ngulo Velocidade angular
0.6 1
0.4 0
dtheta(rad/s)
theta(rad) 0.2 -1
0 -2
-0.2 -3
-0.4 -4
0 1 2 3 4 0 1 2 3 4
t(s) t(s)
posiç ã o Velocidade
1 3
2
0.5
dx(m/s)
1
x(m)
0
0
-1
-0.5 -2
0 1 2 3 4 0 1 2 3 4
t(s) t(s)
53
Simulação – Modelo Linear
%-------------- Simulação do modelo Linear a uma dada condição inicial ------------------
%----------------------------------------- Arquivo 1 -----------------------------------------------
function dx=teste52(t,x)
%Esta função representa matematicamente o pêndulo invertido
%(com linearização) envolvendo o sinal de controle u, sendo:
%x(1)=theta x(2)=theta' x(3)=x x(4)=x' u=-K1x(1)-K2x(2)-K3x(3)-K4x(4) .m
M=2; m=0.1; l=0.5; g=9.8
K1=-131.8861; K2=-28.5152; K3=-30.6122; K4=-23.0892;
dx=zeros(4,1)
dx(1)=x(2)
dx(2)=(M+m)/M/l*g*x(1)+1/M/l*(K1*x(1)+K2*x(2)+K3*x(3)+K4*x(4))
dx(3)=x(4)
dx(4)=-m/M*g*x(1)-1/M*(K1*x(1)+K2*x(2)+K3*x(3)+K4*x(4))
%----------------------------------------- Arquivo 2 ------------------------------------------------
tempo=[0 4] % ou tempo=0:0.001:4; (para fixar um passo)
cond_inicial=[0.6 0 0 0] %cond_inicial=[theta theta' x x']
[T,X]=ode45('teste52',tempo,cond_inicial)
subplot(2,2,1); plot(T,X(:,1))
xlabel('t(s)'); ylabel('theta(rad)'); title ('ângulo')
.m
subplot(2,2,2); plot(T,X(:,2))
xlabel('t(s)'); ylabel('dtheta(rad/s)'); title ('Velocidade angular')
subplot(2,2,3); plot(T,X(:,3))
xlabel('t(s)'); ylabel('x(m)'); title ('posição')
subplot(2,2,4); plot(T,X(:,4))
xlabel('t(s)'); ylabel('dx(m/s)'); title ('Velocidade')
54
Simulação – Modelo Linear
â ngulo Velocidade angular
0.6 1
0.4 0
dtheta(rad/s)
theta(rad)
0.2 -1
0 -2
-0.2 -3
-0.4 -4
0 1 2 3 4 0 1 2 3 4
t(s) t(s)
posiç ã o Velocidade
0.8 3
0.6
2
dx(m/s)
0.4
x(m)
1
0.2
0
0
-0.2 -1
0 1 2 3 4 0 1 2 3 4
t(s) t(s)
55
Comparação entre os Modelos .m
ângulo Velocidade angular
0.5 1
Não-Linear
dtheta(rad/s)
theta(rad)
Linear 0
0 Não-Linear
-1
Linear
-0.5 -2
0 1 2 3 4 0 1 2 3 4
t(s) t(s)
posição Velocidade
0.5 2
Não-Linear Não-Linear
Linear 1 Linear
dx(m/s)
x(m)
0
0
-0.5 -1
0 1 2 3 4 0 1 2 3 4
t(s) t(s)
Comparação entre os modelos para a condição inicial x0=[0.3 0 0 0]
56
Comparação entre os Modelos
ângulo Velocidade angular
1 4
Não-Linear Não-Linear
2
Linear Linear
0.5
0
dtheta(rad/s)
theta(rad)
0 -2
-4
-0.5
-6
-1 -8
0 1 2 3 4 0 1 2 3 4
t(s) t(s)
posição Velocidade
2 6
Não-Linear Não-Linear
1.5 4
Linear Linear
2
1
dx(m/s)
x(m)
0
0.5
-2
0 -4
-0.5 -6
0 1 2 3 4 0 1 2 3 4
t(s) t(s)
Comparação entre os modelos para a condição inicial x0=[0.9 0 0 0]
57
Exemplo – Plano de fase
Obter o plano de fase do sistema representado pela equação
x xx x 0,
para a condição inicial x 3, x 0.
58
Exemplo – Plano de fase
%--------------------------------------- Plano de Fase ----------------------------------------
%----------------------------------------- Arquivo 1 -------------------------------------------
function dx=plano(t,x)
%Plano de Fase da equação: d2x/dt2+[Link]/dt+x=0
dx=zeros(2,1);
dx(1)=x(2); .m
dx(2)=-x(1)*x(2)-x(1);
%----------------------------------------- Arquivo 2 -------------------------------------------
clear;
tempo=0:0.001:9;
cond_inicial=[0 3]; %cond_inicial=[x x']
[T,X]=ode45('plano',tempo,cond_inicial); .m
figure(1)
plot(X(:,1),X(:,2))
xlabel('x')
ylabel('dx/dt')
title ('Plano de Fase')
grid
figure(2);comet(X(:,1),X(:,2)) % para verificar o sentido
59
Exemplo – Plano de fase
Plano de Fase
3
2.5
1.5
dx/dt
0.5
-0.5
-1
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
x
60
3 – Construção do Lugar das
Raízes com MATLAB
61
Sistema de Controle
Lugar das raízes (Root Locus)
62
Exemplo - comando rlocus e rlocfind :
Pólos em Malha Aberta: -4, 0, -2+4i, -2-4i
Zeros em Malha Aberta: -8
>>num=[1 8];
>> den=[1 8 36 80 0];
>> rlocus(num,den)
>> [K,p]=rlocfind(num,den)
63
Exemplo - comando rltool:
Pólos em Malha Aberta: -4, 0, -2+4i, -2-4i
Zeros em Malha Aberta: -8
>> num=[1 8];
>> den=[1 8 36 80 0];
>>[A,B,C,D]=tf2ss(num,den);
>> sys=ss(A,B,C,D);
>> rltool(sys)
.m
64