Carlos André Vaz Junior
[email protected]
http://www.eq.ufrj.br/links/h2cin/carlosandre
EQ/UFRJ
O mundo MATLAB
Mais de 1 milhão de resultados
EQ/UFRJ
Ajuda
?
http://newsreader.mathworks.com
EQ/UFRJ
Livros
EQ/UFRJ
Exemplo 1
Exemplo 1
EQ/UFRJ
Exemplo 1
No processo de compra de um equipamento de
pesquisa, o preço final é a soma do custo do
equipamento (A), frete (B) e impostos (C).
Elaborar uma function que execute o somatório e
retorne o resultado.
EQ/UFRJ
Exemplo 1
EQ/UFRJ
Exemplo 2
Exemplo 2
EQ/UFRJ
Exemplo 2
Tenho uma matriz 2x6 e devo encontrar o elemento
de menor valor. Como se faz?
EQ/UFRJ
Exemplo 2
Achando a posição do menor valor de uma matriz:
x=[1 2 3 4 5 6 ;
2 10 3 3 2 25];
%Forma linear:
xmin=min(x);
xmin=min(xmin);
[i,j]=find(x==xmin);
%Forma condensada:
[i,j]=find(x==(min(min(x))));
EQ/UFRJ
Exemplo 3
Exemplo 3
EQ/UFRJ
Exemplo 3
Para que possamos continuar um projeto é
necessário analisar o comportamento do
parâmetro Z diante das variáveis X e Y.
A equação é:
Z=exp(-0.5*(X2 + Y2))
Elabore um gráfico de Z(x,y).
EQ/UFRJ
Exemplo 3
xx=0:0.01:1;
yy=0:0.01:1;
[X,Y]=meshgrid(xx,yy);
Z=exp(-0.5*(X.^2+Y.^2));
colormap jet
figure(1);surf(X,Y,Z); rotate3d on; shading interp;
EQ/UFRJ
Exemplo 3
EQ/UFRJ
Exemplo 4
Exemplo 4
EQ/UFRJ
Exemplo 4
Nesse exemplo desejamos elaborar um gráfico de superfície
que uma propriedade termodinâmica. O valor dessa propriedade
é função da composição (X1, X2 e K) da mistura. Assim, a base do
gráfico será X e Y, e o eixo vertical representa o valor da
propriedade dado por:
Z=(X.*log(X))+(Y.*log(Y))+(K.*log(K));
Restrição: K refere-se a concentração
do terceiro componente. Desse modo,
K= 1 – X – Y
E K não pode ser negativo!
EQ/UFRJ
Exemplo 4
%Malha triangular da base
%malha da base
xx=0:0.01:1;
yy=0:0.01:1;
[X,Y]=meshgrid(xx,yy);
K=1-X-Y;
%aplica a restrição para usar só a base do triangulo
%onde existe consistência física (o que nao tem vira "Not a Number")
iz=find(K<0);K(iz)=nan;
Z=(X.*log(X))+(Y.*log(Y))+(K.*log(K));
%gráfico da superfície
colormap jet
figure(1);surf(X,Y,Z); rotate3d on; shading interp;
xlabel('X1');ylabel('X2');zlabel('DeltaGi/RT');
EQ/UFRJ
Exemplo 4
EQ/UFRJ
Exemplo 5
Exemplo 5a
EQ/UFRJ
Exemplo 5
EQ/UFRJ
Exemplo 5
Modelagem & Dinâmica de Processos
Modelos simples - o tanque de nível
pode-se escrever o balanço de massa do sistema
Ainda,
dmt 
  FE  F 
dt
(2)
dmt 
dht 
 A
dt
dt
(3)
dht  1
 FE  F 
dt
A
(4)
e, portanto,
1
EQ/UFRJ
Exemplo 5
Modelagem & Dinâmica de Processos
Modelos simples - o tanque de nível
Freqüentemente, considera-se a vazão de saída do tanque
proporcional
à
altura
da
coluna
de
líquido
é inversamente
proporcional a uma resistência ao escoamento (R):
F 
h
R
(5)
Logo,
dht  1 
h
  FE  
dt
A
R
(6)
1
EQ/UFRJ
Exemplo 5
Modelagem & Dinâmica de Processos
Modelos simples - o tanque de nível
Este modelo simples de um tanque de nível, sem balanço de
energia, possui uma solução analítica:
t



ht   RFE 1  e RA 


(7)
Para simular este modelo, basta escolher os valores das
constantes R, A e FE, das condições iniciais h0 e t0.
A simulação da solução analítica do modelo do tanque de
nível é mostrada a seguir.
1
EQ/UFRJ
Exemplo 5
% Definição das constantes do modelo
R = 1;
% h/m2
A = 2;
% m2
Fe = 10;
% m3/h
% Tempo de simulação
t = 0.0 : 0.01 : 10.0;
%h
% Simulação da altura de líquido
h = R*Fe*(1 - exp(-t/(R*A)));
%m
% Visualização da simulação
plot(t,h);
title('Simulação do tanque de nível');
xlabel('Tempo (h)');
ylabel('Altura (m)');
EQ/UFRJ
Exemplo 5
Verifique a consistência do calculo: a
matriz “h” gerada também deve ser 1x1000, já que
cada instante “t” gerou um valor “h”. É sempre útil
conferir a dimensão das variáveis, principalmente a
medida que as rotinas forem tornando-se
complexas.
Dica!
EQ/UFRJ
Exemplo 5
Exemplo 5b
EQ/UFRJ
Exemplo 5
Modelagem & Dinâmica de Processos
Modelos
- o tanque
de nível
Muitas
vezessimples
é muito trabalhoso,
ou mesmo
impossível,
encontrar a solução analítica para o conjunto de
Freqüentemente, considera-se a vazão de saída do tanque
equações
diferenciais. Nesse caso temos que simular
proporcional à altura da coluna de líquido é inversamente
usando solução numérica das equações diferenciais.
proporcional
a uma
resistência
ao escoamento
(R): 1 não tivesse
Vamos
assumir
que
o modelo
do exemplo
solução analítica, e então usar o Matlab para estudar o
h
 do nível com o tempo. A (5)
comportamento da Faltura
R
equação diferencial
será:
Logo,
dht  1 
h
  FE  
dt
A
R
(6)
1
EQ/UFRJ
Exemplo 5
% Definição das constantes do modelo
R = 1;
% h/m2
A = 2;
% m2
Fe = 10; % m3/h
% Tempo de simulação
t = 0.0 : 0.01 : 10.0; % h
% Simulação da altura de líquido
[t,h] = ode45('dhdt',t, 0,[],[R A Fe]);
% Visualização da simulação
plot(t,h);
title('Simulação do tanque de nível');
xlabel('Tempo (h)');
ylabel('Altura (m)');
EQ/UFRJ
function dh = dhdt(t,h,flag,par)
R = par(1);
A = par(2);
Fe = par(3);
dh = (Fe-(h/R))/A;
Exemplo 5
Nesse caso temos uma equação diferencial, então deveremos
usar uma função Matlab específica para a resolução de eq.
diferenciais. No caso temos a ODE45. A função ODE45
implementa um esquema de solução de sistemas de EDO’s por
método de Runge-Kutta de ordem média (consulte o help sobre
ODE45 para maiores detalhes).
[t,h] = ode45('dhdt',t, 0,[],[R A Fe]);
EQ/UFRJ
Exemplo 5
Os parâmetros enviados entre parênteses são aqueles que
devemos passar para a ODE45:
-1º argumento de ode45 é uma string contendo o nome do arquivo .m
com as equações diferenciais. Neste caso, o arquivo chama-se dhdt.m.
-2º argumento é um vetor que pode conter (i) dois elementos: os
tempos inicial e final da integração, ou (ii) todos os valores de tempo para
os quais deseja-se conhecer o valor da variável integrada.
-3º argumento é o vetor contendo as condições iniciais das variáveis
dependentes das EDO’s. Os valores dos elementos do vetor de condições
iniciais precisam estar na mesma ordem em que as variáveis
correspondentes são calculadas na função passada como 1º argumento
para ode45 (neste caso, dhdt.m). Nesse caso em particular só temos uma
variável dependente, assim temos uma única condição inicial.
EQ/UFRJ
Exemplo 5
-4º argumento é o vetor de opções de ode45. Há várias opções do
método que podem ser ajustadas. Entretanto, não deseja-se alterar
os valores-padrão. Neste caso, é passado um vetor vazio, apenas para
marcar o lugar das opções.
-5º argumento é um vetor contendo parâmetros de entrada para a
função dhdt.m. Observe que a função .m deve ler esses parâmetros
na ordem correta (recebe como variável local “par”).
Os resultados da simulação são obtidos nos dois parâmetros
entre colchetes (t , h).
EQ/UFRJ
Exemplo 5
A codificação do arquivo .m segue o mesmo formato já explicado para
funções porém com algumas particularidades.
No caso específico de um arquivo .m que deve ser chamado por uma
função de solução EDO’s (todas as ODExx), a declaração deste arquivo
deve seguir a sintaxe:
function dy = nomefun(t, y, flag, arg1, ..., argN)
onde
•dy é o valor da(s) derivada(s) retornadas
•t e y são as variáveis independente e dependente, respectivamente.
•Opcional: caso deseje-se receber outros parâmetros, a função deve
receber um argumento marcador de lugar chamado flag. Após este,
ela recebe quaisquer outros parâmetros.
EQ/UFRJ
Exemplo 5
Exemplo 5c
EQ/UFRJ
Exemplo 5
É necessário agora guardar alguns termos intermediários
realizados no cálculo da function. Por exemplo:
function dh = dhdt(t,h,flag,par)
global Y tvetor
R = par(1);
A = par(2);
Fe = par(3);
dh = (Fe-(h/R))/A;
Quero guardar:
Y=Fe+A+dh;
Como se faz?
EQ/UFRJ
Exemplo 5
O programa principal
ganhou apenas três linhas
novas. O comando global
vai estabelecer comunicação
com a function e os comandos
figure e plot fazem o gráfico
do novo parâmetro.
EQ/UFRJ
% Definição das constantes do modelo
clear all
global Y tvetor
R = 1;
% h/m2
A = 2;
% m2
Fe = 10; % m3/h
% Tempo de simulação
t = 0.0 : 0.01 : 10.0; % h
% Simulação da altura de líquido
[t,h] = ode45('dhdt',t, 0,[],[R A Fe]);
% Visualização da simulação
plot(t,h);
title('Simulação do tanque de nível');
xlabel('Tempo (h)');
ylabel('Altura (m)');
figure(2)
plot(tvetor,Y,'*b')
Exemplo 5
Dentro da function é necessário
usar o comando global para estabelecer
a troca dos novos dados.
O novo dado é: Ylinha=Fe+A+dh;
Para acumular usamos a matriz Y.
O vetor “tvetor” guarda o tempo
referente para cada valor de Y
calculado.
EQ/UFRJ
function dh = dhdt(t,h,flag,par)
global Y tvetor
R = par(1);
A = par(2);
Fe = par(3);
dh = (Fe-(h/R))/A;
Ylinha=Fe+A+dh;
Y=[Y Ylinha];
tvetor=[tvetor t];
Exemplo 5
EQ/UFRJ
Exemplo 6
Exemplo 6
EQ/UFRJ
Exemplo 6
EQ/UFRJ
Exemplo 6
Modelagem & Dinâmica de Processos
Modelagem & Dinâmica de Processos
Modelos simples - tanque de aquecimento
Modelos simples - tanque de aquecimento
Como no caso anterior, o balanço de massa pode ser escrito
   V dT
d :VT
como
dt
dt
T
dV
dh 
 dT
 A h
T

dt
dt 
 h dt
1
dht 
  FE  
dt
A
R
 dT T 
h 
C p A  h
  FE     FE HE  FH   Q
dt
A
R

 como:
O
balanço de energia
é escrito
 UTFh H FH
dT C1 d FVT
FE   Q
U 
E TE  
 p


E 
E T
 A

dt
h  dtA
C p 

C
p 


(9)
(6)
(10)
(8)
(11)
1
EQ/UFRJ
1
Exemplo 6
Traduzindo as equações diferenciais para o Matlab:
EQ/UFRJ
Matlab
Real
dy(1)
dh/dt
y(1)
h
dy(2)
dT/dt
y(2)
T
Exemplo 6
% Definição das constantes do modelo
R = 1;
% h/m2
A = 2;
% m2
Fe = 10;
% m3/h
Cp = 0.75; % kJ/(kg . K)
Ro = 1000; % kg/m3
U = 150;
% kJ/(m2 . s . K)
Te = 530;
%K
Th = 540;
%K
% Tempo de simulação
t = 0.0 : 0.01 : 10.0; % h
% Simulação do modelo
[t,y]=ode45('dydt',t,[(5/A) Th],[],[U A Ro Cp Fe R Te Th]);
EQ/UFRJ
Exemplo 6
% Visualização da simulação
figure(1);
plot(t,y(:,1));
title('Tanque de aquecimento');
xlabel('Tempo (h)');
ylabel('Altura (m)');
figure(2);
plot(t,y(:,2));
title('Tanque de aquecimento');
xlabel('Tempo (h)');
ylabel('Temperatura (K)');
EQ/UFRJ
Exemplo 6
A única modificação em relação ao exemplo anterior é que
estamos passando duas condições iniciais (pois existem duas
variáveis dependentes):
[t,y]=ode45('dydt',t,[(5/A) Th],[],[U A Ro Cp Fe R Te Th]);
EQ/UFRJ
Exemplo 6
A função .m tem o código apresentado a seguir:
function dy = dydt(t,y,flag,par);
U = par(1);
A = par(2);
Ro = par(3);
Cp = par(4);
Fe = par(5);
R = par(6);
Te = par(7);
Th = par(8);
dy(1) = (Fe-(y(1)/R))/A;
dy(2) = (1/y(1))* ( ((Fe*Te/A)+(U*Th/(Ro*Cp)))...
- ( y(2)*((Fe/A)+(U/(Ro*Cp)))) );
dy = dy(:);
EQ/UFRJ
Exemplo 6
O vetor dy é criado como vetor linha (dy(1)) e
(dy(2)). Porém temos que retornar como vetor coluna.
Use o comando:
matriz coluna = matriz linha (:)
Dica!
EQ/UFRJ
Exemplo 6
Quando for fazer os gráficos no programa principal
lembre-se que a primeira coluna de “dy” refere-se a “h” e a
segunda a “T”. Então para graficar h vs. tempo faça:
figure(1);
plot(t,y(:,1));
title('Tanque de aquecimento');
xlabel('Tempo (h)');
ylabel('Altura (m)');
Dica!
EQ/UFRJ
Exemplo 7
Exemplo 7
EQ/UFRJ
Exemplo 7
Em uma amostragem experimental obtivemos o seguinte
conjunto de pontos:
xexp=[0 1 2 3 4 5];
yexp=[-25 -10 7 6.7 -5 -15];
Ajuste um polinômio de 2ª ordem para os dados
EQ/UFRJ
Exemplo 7
Programa principal:
close all
% chute inicial:
A=4;
B=2;
C=3;
X0=[A;B;C];
options = optimset('display','iter');
X = fminsearch('UUm',X0,options);
disp('Resultado (A, B, C)')
disp(X)
figure(2)
xexp=[0 1 2 3 4 5];
yexp=[-25 -10 7 6.7 -5 -15];
x=0:0.01:5;
y=(X(1)*(x.^2))+(X(2).*x)+X(3);
plot(x,y,'-b')
hold on
plot(xexp,yexp,'*r')
EQ/UFRJ
Exemplo 7
Função UUm:
function [erro] = UUm(entra)
A=entra(1);
B=entra(2);
C=entra(3);
xexp=[0 1 2 3 4 5];
yexp=[-25 -10 7 6.7 -5 -15];
ymod=(A*(xexp.^2))+(B.*xexp)+C;
erro=(yexp-ymod).^2;
erro=sum(erro);
plot(xexp,yexp,'*r')
hold on
plot(xexp,ymod,'-b')
pause(0.1)
EQ/UFRJ
Exemplo 8
Exemplo 8
EQ/UFRJ
Exemplo 8
A seguinte malha de controle foi elaborada no Simulink.
Usar o Matlab para ajustar o controlador.
degrau unitário
no instante 5
EQ/UFRJ
P
I
D
Exemplo 8
Programa principal:
clear all
close all
warning off
options = optimset('display','iter');
global P I D erro
Pmin = fminsearch('custo', [1 5 1],options)
EQ/UFRJ
Exemplo 8
Função “custo”:
function [erro] = custo(x)
global P I D erro
P=x(1);
I=x(2);
D=x(3);
[T]=sim('malha',[0 65]);
erro=sum(erro.^2);
EQ/UFRJ
Exemplo 8
Solução encontrada para degrau unitário no SP:
Pmin = 4.5075 2.6329 -0.0000
EQ/UFRJ
Exemplo 8
Arquitetura Simulink usada para gerar o gráfico
do slide anterior:
EQ/UFRJ
Exemplo 8b
Exemplo 8b
EQ/UFRJ
Exemplo 8b
Ao invés de minimizar o somatório quadrático do erro,
posso minimizar o somatório quadrático ponderado com
o tempo. Ou seja, erros em tempos mais elevados são mais
significativos.
EQ/UFRJ
Exemplo 8b
EQ/UFRJ
Exemplo 8b
Programa principal:
clear all
close all
warning off
options = optimset('display','iter');
global P I D erro tempo
Pmin = fminsearch('custo', [10 5 1],options)
EQ/UFRJ
Exemplo 8b
Função custo:
function [erro] = custo(x)
global P I D erro tempo
P=x(1);
I=x(2);
D=x(3);
[T]=sim('malha',[0 65]);
%erro=sum(erro.^2); % somatorio quadratico do erro
erro=sum((erro.*tempo).^2); % somatorio quadratico
do erro ponderado com o tempo
EQ/UFRJ
Exemplo 8b
Resultado obtido:
Pmin = 25.8333
5.3333 -0.0000
Chute inicial usado:
10 5 1
EQ/UFRJ
Exemplo 9
Exemplo 9
EQ/UFRJ
Exemplo 9
Um sistema de detecção de falhas retorna em seu relatório
o índice de confiança de cada um dos sensores disponíveis.
O resultado está na forma de texto:
SVI-D sensor 1:0.87514
SVI-D sensor 2:0.96246
SVI-D sensor 3:0.99978
SVI-D sensor 4:0.942
SVI-D sensor 5:0.86938
SVI-D sensor 6:0.98819
SVI-D sensor 7:0.94546
SVI-D sensor 8:0.82448
SVI-D sensor 9:0.99934
SVI-D sensor 10:0.91721
SVI-D sensor 11:0.99358
SVI-D sensor 12:0.96833
SVI-D sensor 13:0.92814
SVI-D sensor 14:0.83297
SVI-D sensor 15:0.97261
SVI-D sensor 16:0.98094
EQ/UFRJ
SVI: índice de validade do sensor
Exemplo 9
A questão então é achar o sensor que apresenta o menor
índice de validade. Ou seja, achar o menor SVI. Como fazer
isso sem precisar digitar tudo de novo?
EQ/UFRJ
Exemplo 9
celula=importdata('texto.m');
texto=char(celula);
resultado=[];
sensores=[];
for i=1:16,
posicao=find(texto(i,:)==':');
valor=texto(i,posicao+1:end);
valor=str2num(valor);
resultado=[resultado;valor];
end
for i=1:16,
sensor=find(ordenado(i,:)==resultado);
% Ordena do menor para o maior:
ordenado=sort(resultado);
sensores=[sensores;sensor];
end
disp('Valores ordenados')
ordenado
disp(' ')
disp('Em ordem de sensores:')
sensores
EQ/UFRJ
Exemplo 10
Exemplo 10
EQ/UFRJ
EQ/UFRJ
EQ/UFRJ
EQ/UFRJ
EQ/UFRJ
EQ/UFRJ
EQ/UFRJ
EQ/UFRJ
EQ/UFRJ
Exemplo 11
Exemplo 11
EQ/UFRJ
Exemplo 11
Na compra de uma calculadora gráfica, a loja ofereceu duas propostas
de financiamento – proposta “A” e “B”. A proposta A é composta
por 7 parcelas mensais iguais de 114 reais cada. Já a proposta B prevê
10 parcelas mensais iguais de 98 reais cada. Qual é a melhor opção de
compra considerando a taxa de juros oferecida em investimentos
denominados “renda fixa”?
A princípio poderia resolver o problema simplesmente multiplicado
114 x 7 e 10 x 98, achando o valor final pago. Os valores
encontrados seriam 798 e 980. Logo, a Proposta A parece mais
favorável para o comprador.
EQ/UFRJ
Exemplo 11
É importante lembrar, porém, que essa forma de resolução não
considera que o dinheiro desvaloriza-se ao longo dos meses. Ou seja,
o poder de compra de 100 reais hoje, é superior ao poder de compra
de 100 reais daqui a 10 meses. Outra forma de pensar é considerar o
“custo de oportunidade” – a taxa de retorno livre de que conseguiria
para o meu dinheiro caso, ao invés de pagar agora, investisse. De
uma forma ou de outra o que precisamos é do VALOR PRESENTE
(VP) de cada série de pagamentos, sendo os pagamentos descontados
a dada taxa de juros. Para trazer VALOR FUTURO (VF) para valor
presente usa-se a fórmula:
EQ/UFRJ
Exemplo 11
Onde “i” é a taxa de juros mensal e “n” o número de meses
entre o VF e o VP.
O código MATLAB abaixo representa os valores presentes para
cada taxa de juros implementada. O programa aponta ainda o
valor de juros que o VP de ambas as séries tornam-se iguais.
Depois, apresenta-se detalhadamente a explicação do código
implementado.
EQ/UFRJ
Exemplo 11
clc
close all
clear all
ivetor=0:0.01:0.50;
VPvetor114=[];
VPvetor98=[];
prompt{1}='Número de meses do pagamento da serie A:';
prompt{2}='Número de meses do pagamento da serie B:';
prompt{3}='Valor de cada parcela da serie A:';
prompt{4}='Valor de cada parcela da serie B:';
resposta=inputdlg(prompt,'Calculo da taxa de equilibrio');
nummeses114=str2num(char(resposta(1)));
nummeses98=str2num(char(resposta(2)));
v114=str2num(char(resposta(3)));
v98=str2num(char(resposta(4)));
EQ/UFRJ
Exemplo 11
for J = 1:length(ivetor),
i=ivetor(J);
VP=[];
for K = 1:nummeses114,
VP(K)=v114/(1+i)^K;
end
VPfinal=sum(VP);
VPvetor114=[VPvetor114, VPfinal];
VP=[];
for K = 1:nummeses98,
VP(K)=v98/(1+i)^K;
end
VPfinal=sum(VP);
VPvetor98=[VPvetor98, VPfinal];
end
EQ/UFRJ
Exemplo 11
plot(ivetor*100,VPvetor114,'-b')
hold on
plot(ivetor*100,VPvetor98,'-r')
title('Valor presente das parcelas a serem pagas')
legend( [ num2str(nummeses114), ' parc de ', num2str(v114) ,' reais cada'] , ...
[ num2str(nummeses98), ' parc de ', num2str(v98) ,' reais cada' ] )
xlabel('Taxa de juros mensal')
ylabel('Valor presente em Reais')
if (VPvetor114(1)<VPvetor98(1)),
posicoes=VPvetor114<VPvetor98;
achaZero=find(posicoes==0);
achaPrimeiroZero=min(achaZero);
plot(100*ivetor(achaPrimeiroZero),VPvetor114(achaPrimeiroZero),'ok')
else
posicoes=VPvetor98<VPvetor114;
achaZero=find(posicoes==0);
achaPrimeiroZero=min(achaZero);
plot(100*ivetor(achaPrimeiroZero),VPvetor114(achaPrimeiroZero),'ok')
end
EQ/UFRJ
text(100*ivetor(achaPrimeiroZero),50+VPvetor114(achaPrimeiroZero), ...
['Juros de equilibrio (a.m.) = ',num2str(100*ivetor(achaPrimeiroZero)),' %'] )
Exemplo 11
EQ/UFRJ
Carlos André Vaz Junior
[email protected]
http://www.eq.ufrj.br/links/h2cin/carlosandre
EQ/UFRJ