107506– Laboratório de Controle de Processos

Prof. Eduardo Stockler Tognetti

Departamento de Engenharia Elétrica

Universidade de Brasília

2º Semestre 2017

Aula: Matlab

Introdução

• Ao contrário de linguagens clássicas como C e Fortran, no ambiente MATLAB o usuário não se preocupa com itens como declaração de variáveis, alocação de memória, utilização de ponteiros e outras tarefas de rotina.

• O MATLAB apresenta uma série de funções matemáticas já implementadas que podem ser utilizadas em uma rotina construída pelo usuário. Estas funções são agrupadas de acordo com a área de interesse em toolboxes e armazenadas em diretórios específicos. Qualquer função a ser utilizada deve estar no diretório de trabalho ou caminho (path) do MATLAB.

Ambiente

Editor de códigos (scripts), arquivo .m salvo

em disco

Caminho (path) do arquivo .m

Janela de comando, executa

diretamente os comandos, não salva em disco

Lista de variáveis na memória da

área de trabalho (workspace)

Lista os arquivos do diretório do

“path”

Histórico dos comandos

executados no Command Window

Inicia Simulink

Comandos Básicos Comando Descrição Exemplo

help ajuda help plot

who, whos informa as variáveis presentes no workspace whos

which informa o diretório de execução de uma determinada função

which ode45

lookfor Procura uma palavra na biblioteca de funções lookfor exponential

load, save carrega / salva em disco as variáveis do workspace como arquivo de dados (extensão .mat)

save dados.mat

clc, close Limpa a tela, figuras clc; close all

clear limpa variáveis do workspace clear a; clear all;

diary salva o texto de uma seção do MATLAB no diretório corrente

diary texto.txt

figure cria uma figura nova figure; figure newfig

cd Muda ou informa o diretório corrente cd ..; cd c:\temp; cd

Operações e Variáveis

• Operações aritméticas padrão ou pré-definidas:

» (5+4)/(6*2)

» sin(pi/4)^2 – 1

• Nome de variáveis deve começar com letra:

» raio3 = 2;

» circunferencia = 2*pi*raio3;

Obs.: “;” no final omite impressão na tela.

Funções Elementares

Vetores e matrizes

Atribuindo valores:

• v = [1 2 3 4]; v=1:4; v=[1;2;3;4]’; A = [1 2; 3 4];

• x=(0:0.1:1)*pi; y= linspace(0,pi,11); z= logspace(0,2,11);

• I=eye(3); U=ones(3); U=ones(3,2); Z=zeros(3); Z=zeros(3,1); A=rand(2,4);

Atribuindo elementos específicos:

• v(2)=20; A(2,1)=30; v(1:end)=1;

• A(:,1)=pi; A(2,:)=log(exp(pi)); A([1 3], : )=sin(pi);

Informações:

• [l,c] = size(A); c = length(v); display(A);

Vetores e matrizes

Operações: • (+, -, *, /, ^ ) : operações matriciais • (.*, ./, .^ ) : operações elemento a elemento • eig(A), norm(A), det(A), inv(A), rank(A), null(A), expm(A) auto valores/vetores, norma, determinante, inversa, posto da matriz A, espaço nulo e exp. de matriz, respectivamente. • diag(A), diag([1 2]): extrai a diagonal da matriz A, gera uma

matriz diagonal com os elementos 1 e 2. • x=A\b solução de A*x=b (≡ inv(A)*b, se A quad. ñ sg.) • x=b/A solução de x*A=b (≡ b*inv(A), se A quad. ñ sg.)

help matfun: help on matrix functions (numerical linear algebra)

Vetores e matrizes

Matrizes multidimensionais: • M(1,1,1,1)=1; M=randn(2,2,3); • Obs.: Converter matriz 3D para 1D: vet = squeeze(M); Células: Criando: • a=1:3; b=[5 6;7 8]; c={a b}; d.a=a; d.b=b; • z={a b c d}; • C={[1 2],[4 5;6 7]}; Mostrando ou atribuindo valores: • z{2}=sin(pi); z{2}(1,2); z{3}{2}(1,2); z{3}{2}(2,2)=80; celldisp(z) Estruturas: • dados.nome='prog'; dados.mat={[1] [1;2]}; • s = struct('Disciplina',lab','Nota',5.0); • valores = estrutura.signal.values;

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);

Exemplo 2

Exemplo: Achando a posição do menor valor de uma matriz:

%Forma condensada: [i,j]=find(x==(min(min(x))));

Vetores e matrizes

%Criando células (cells) temperature(1,:) = {'2009-12-31', [45, 49, 0]}; temperature(2,:) = {'2010-04-03', [54, 68, 21]}; temperature(3,:) = {'2010-06-20', [72, 85, 53]}; allTemps = cell2mat(temperature(:,2)); dates = datenum(temperature(:,1)); plot(dates,allTemps)

Exemplo 2

Exemplo: Manipulando células (cell)

%Criando células: for i =1:10 c{i} = i; end

Vetores e matrizes

%Célula para matriz: C = {[1] [2 3 4]; [5; 9] [6 7 8; 10 11 12]} celldisp(C) M = cell2mat(C)

%Criando células: c = cell(1,3); c{1} = 1; c{2} = [1 2 3]; c{3} = [1 2; 3 4];

Atribuindo valores e imprimindo Atribuindo múltiplos valores:

• [A,B,C,...] = deal(X,Y,Z,...); [A,B,C,...] = deal(X{:}); • [A,B,C,...] = deal(S.FIELD);

Exemplo:

• sys = {rand(3) ones(3,1) eye(3) zeros(3,1)}; %cell

• [a,b,c,d] = deal(sys{:})

• z = input ( ‘Valores de z: ’); • temp=78; disp(temp); disp (‘graus Celsius’); • display(['A temperatura é ' num2str(temp) ' graus Celsius']) • fprintf (formato, matriz);

– fprintf (‘A temp. é %4.1f graus Celsius \n’, temp)

Strings

Entrada: • s = 'Este texto é uma string'; display(s(1:4)); Manipulação: • eval(S): avalia a string S como uma expressão • str2num('7'): converte a string 7 no número 7 • num2str(7): converte o número 7 em string • strcat(S,T): concatena as strings S e T • strcmp(S,T): compara as strings S e T Exemplos: • eval('R=1, A=2, h0=10'); save(['Nome_' date]); • nome='a'; eval(['P',nome,' = 10;']); P = eval(['P',nome]) help strfun: help on character strings.

Polinômios

• Polinômio:

• Exemplo:

• Operações: polinômios com coefs. x, y e raízes r – Retorna as raízes: roots(x)

– Retorna os coeficientes: poly(r)

– Avaliar polinômio x em s (ex.: s=5): polyval(x,s)

– Multiplicação de polinômios: z=conv(x,y)

– Divizão de polinômios: [z,t]=deconv(x,y)

– Derivar polinômios: polyder(x)

– Soma de polinômios: x+y (completar c/ zeros à esquerda)

Polinômios

• Exemplo: Ajuste polinomial por polinômio de 2ª ordem x=[0 1 2 3 4 5]; y=[4 5 7 6 6 5]; %acha os parametros do polinomio [p,s]=polyfit (x,y,2); %plota o resultado xx=0:0.01:5; pop4 = polyval (p,xx); plot (x,y,'b.') hold on plot (xx,pop4,'r -')

Integração e Derivação Numérica

Dada uma função utilizar integral. Ex.: Q = integral(@(x) exp(-x.^2).*log(x).^2,0,Inf)

Dado um vetor de dados utilizar trapz. Ex.: Q = trapz([1 4 9 16 25])

Integração e Derivação Numérica

Dada uma função utilizar “integral”. Ex.: Q = integral(@(x) exp(-x.^2).*log(x).^2,0,Inf)

Dado um vetor de dados utilizar “trapz”. Ex.: Q = trapz([1 4 9 16 25])

• Obs.: “cumtrapz” fornece o valor acumulado.

Integração numérica:

h = 0.001; % step size

X = -pi:h:pi; % domain

f = sin(X); % range

Y = diff(f)/h; % first derivative

Z = diff(Y)/h; % second derivative

plot(X(:,1:length(Y)),Y,'r',X,f,'b', X(:,1:length(Z)),Z,'k')

Derivação numérica:

Matemática simbólica Definido variáveis simbólicas » syms x a b Construindo uma função simbólicas com as variáveis criadas » z = 1/(a+b*(cos(x))^2) Mostrando a função » pretty(z); Simplificando » simple(z) Diferenciando em relação a x » dz = diff(z,x) Suvbstituindo variáveis por valores » a=0.5; b=1; » dzab=subs(dz) » z2 = subs (z,'a',100) Resolvendo equações » pts_crits=solve(dzab) Ajuda: • Outros: solve (resolve eqs. algébricas); dsolve (resolve edo’s); fzero; fsolve; compose • Otimização: optimset; fmincon; fminsearch; lsqnonlin

Matemática simbólica

Curve Fitting

• Exemplo: f = @(x) -5*x.^2+exp(-2*x) + 0.1*randn; fplot(f,[-2 2]) %Data x = linspace(-2,2,1e3)'; for i=1:length(x) y(i,1) = f(x(i)); end type = 'exp2'; %poly3, smoothingspline (sinal ruidoso)... [curve, goodness] = fit( x, y, type); % polyfit() cftool(x,y)

Otimização Fonte: Mathworks

Otimização Graphical optimization tool: Optimization App Sintaxe: optimtool Fonte: https://www.mathworks.com/help/optim/ug/graphical-optimization-tool.html

Gráficos 2D Gráficos 2D: » x = [-pi:0.1:pi] ; y = cos(x); plot(x,y); Duas curvas no mesmo gráfico: » plot(x, 1./x, 'r', x, y, '--b'); ou » figure; hold on; plot(x, 1./x,'r'); plot(x, y,'k'); hold off; grid; Curvas em gráficos diferentes mas na mesma janela: » subplot(1,2,1);plot(x, 1./x,' *b' );subplot(1,2,2);plot(x,y,'--m' )

Propriedades: • title; xlabel; ylabel; legend; grid; ginput; gtext; • axis([XMIN XMAX YMIN YMAX]) • fig01 = figure; saveas(fig01,’nomefig','jpeg'); • print -depsc -tiff -r300 filename help plot

Gráficos 3D

» figure; % 3-D mesh surface » x = linspace(-10,10,100)'; » y = x; » [X,Y] = meshgrid (x,y); » Z = X.^2 + Y.^2; » mesh(Z);

» figure; % Plot lines and points in 3-D space » t = 0:pi/50:10*pi; » plot3(sin(t),cos(t),t);

• Outros (2D): stem, line, semilogx, loglog, scatter, bar; area; pie; hist; • Outros (3D): plot3, stem3, surf, surfl; contour; quiver; fill

Exemplo: duas curvas no mesmo gráfico

x = [-2*pi:0.1:2*pi] ; y = cos(x); y2 = sin(x); %Exemplo1 figure plot(x,y,‘-b',x,y2,'r') %Propriedades xlabel('t') ylabel('y(t)') title('Resposta do processo') legend('y1','y2') grid

%Exemplo2 figure plot(x,y,'b') hold on plot(x,y2,'r') hold off

%Exemplo3 figure subplot(2,1,1) plot(x,y,'b') subplot(2,1,2) plot(x,y2,'r')

Gráfico da derivada numérica de um sinal

h = 0.1; %passo

x = -2*pi:h:2*pi;

y = cos(x);

plot(x,y)

y2 = diff(y)/h; %y2(k) = ( y(k+1) – y(k) )/h

%Gráfico de y e gráfico da derivada de y com ajuste das dimensões dos vetores

plot(x(1:end-1),y(1:end-1),'-b',x(1:end-1),y2,'r')

Exercício 1

2. Seja a célula C = {[1] [2 3 4]; [5; 9] [6 7 8; 10 11 12]}

a) Visualize o conteúdo (comando “celldisp”) b) Modifique o valor do elemento (1,1) da célula de “1” para “10” c) Passe todos os elementos de C para uma matriz M (comando “cell2mat”) d) Passe cada um dos elementos da célula para M1, M2, M3 e M4

(manualmente e depois via comando “deal”) e) Seja N(1,1,1:4)=rand. Converta N para um vetor v (comando “squeeze”)

1.

Exercício 2

1. Plote num mesmo gráfico (diferencie as curvas por cores e marcadores) as três funções abaixo utilizando (i) Comando plot (ii) Comando hold

2. a) Resolva g(x)=0 (ache as raízes x), em que g(x) = 16/x - (x-5)^2 = 0

b) Resolva g2(x) = dg(x)/dx = 0 c) Plote g2(x) para -10 < x < 10 Dica: para obter valores numéricos use o comando “eval”

Leitura e escrita de arquivos

Lendo

• arquivo = fopen(’temperaturas.txt’);

• A = fscanf(arquivo,’Hora %d, Temperatura %f\n’,[2 7])

• A = fscanf(arquivo,’Hora %d, Temperatura %f\n’,[2 Inf])

• fclose(arquivo); disp(A)

Escrevendo

• arquivo = fopen(’info.txt’,’w’)

• fprintf(arquivo,’%d’,1); fclose(arquivo)

Obs.: No Excel: [dados,texto,resto] = xlsread(arquivo,planilha)

Script-files

% Programa Script .m: % Descrição do programa % Inicialização clc, clear all , clear global % Dados do programa perc=[100]; % Entrada de dados entrada=input('Entre com os dados [a,b] '); % Execução Entradaperc = roundn(entrada/perc,-2)

Definindo Funções

Definindo uma função anônima (function handle):

» y = 7;

» fh = @(x)x.^2+y;

» z = fh(2)

ou

» funcao = @(x)cos(x)*x^2;

Plotando

» fplot(funcao,[-2*pi 2*pi]);

Funções em arquivos .m

Formato:

function [output_args] = fct_name(input_args)

[y1, y2,…,yN] = funcao(u1,u2,…,uM)

• Ex: Arquivo “funcao.m”

A função “funcao.m” pode ser chamada a partir de outro arquivo.m ou do Command Windows

Ex.: >> [y1, y2] = funcao(4,7);

function [y1, y2] = funcao(u1,u2) y1 = u1 + u2; y2 = u1*u2 end

Funções em arquivos .m

Formato: • function [output_args] = fct_name(input_args) Exemplo: 1. Definir o diretório no “path” 2. Criar arquivo “media.m” 3. Escrever o código abaixo no arquivo “media.m”

4. Executar na linha de comando (Command Window) » z = 1:99; » average(z)

function y = media(x) if ~isvector(x) error('Input must be a vector') end y = sum(x)/length(x); end

Funções em arquivos .m

• As funções podem estar dentro de outra função

• Ex.: Seja o arquivo “exemplo.m”

Ex.: >> valor = exemplo(4,7);

function y = exemplo(u1,u2) x1 = u1 + u2; x2 = u1*u2; y = media(x1,x2); function saida = media(a,b) saida = (a+b)/2; end end

Funções em arquivos .m

• Funções podem não ter argumentos de entrada

• Nested functions pode usar variáveis que não são explicitamente usadas como argumentos de entrada

Ex.: >> valor = exemplo

function y = exemplo %parent function u = rand(1,2); x1 = u(1) + u(2); x2 = u(1)*u(2); y = media; function saida = media %nested function saida = (x1+x2/2); end end

Function-files

• function [output_args] = fct_name(input_args) function [saida] = funcao(a,b) % function [saida] = funcao(a,b) % Descrição da função % Dados do programa perc=[100]; % Execução saida = roundn([a b]/perc,-2) • Variáveis globais: >> global variavel • Outros comandos: return (functions), break (for or while loops), continue

(for or while loops):

Uso de variáveis globais

function [saida] = nomefuncao(a,b) global perc perc=100; saida = conta([a b],perc); function y = conta(x) global perc y = roundn([x(1) x(2)]/perc,-2); end end Outros comandos: return (functions), break (for or while loops), continue (for or while

loops):

Comandos usados em functions

• return: Immediately exits the function in which it appears. Control passes to the caller of the function

• break (for or while loops): Exits the loop in which it appears. In nested loops, control passes to the next outer loop.

• continue (for or while loops): Skips any remaining statements in the current loop. Control passes to next iteration of the same loop.

Exemplo: variáveis globais

function [y1,y2] = minhafuncao(x1,x2) global a % variável global que pode ser utilizada em outras instâncias do código a = 3; y1 = x1*x2 + a; y2 = x1 + novafuncao(x2); % é possível criar uma função dentro de um função function y = novafuncao(u) global a y = 10*u + a;

Criar arquivo minhafuncao.m com o seguinte código:

Executar função ‘minhafuncao’ em novo arquivo .m ou no Command Window: [saida1,saida2] = minhafuncao(2,3)

Laços

The ‘for’ loop:

for ‘counter’ (statements) end

The ‘while’ loop:

while ‘conditional statement’ (statements) end

The ‘if’ statements: if ‘variable’ ‘conditional statement’ constant (statements) elseif ‘variable’ ‘conditional statement’ constant (statements) else (statements) end

Relações

Operadores relacionais: Operadores lógicos:

Obs.: ‘=’ atribui valor e ‘==’ estabelece uma relação.

a = 5; b = 2;

if a == 5 & b == 2 display(‘Acertou!'); end

Laços for i=1:10 if i <= 5 c(i)=i+1; elseif i < 7 c=[c i]; else c=c.^2; end end method = 'Bilinear'; switch lower(method) case {'linear','bilinear'} disp('Method is linear') case 'cubic' disp('Method is cubic') otherwise disp('Unknown method.') end

test=zeros(3,3); k=0; for i=1:10 while k == 0 test=test+eye(3); if trace(test) < 20 k=0; else k=1; end end end Obs.: ‘=’ atribui valor e ‘==’ é uma relação.

a = 5; if a == 5 display(‘a vale 5'); end

Observação

Método 1: versão com laço

x=0;

for i=1:1000

y(i)=log10(x);

x=x+0.01;

end

Método 2: versão vetorial

x=0:0.01:1000;

y=log10(x);

• Prefira Método 2 ao Método 1.

Exercício 3 Implementar:

1. Define two functions in a file named stat2.m, where the first function calls the second.

function [m,s] = stat2(x) n = length(x); m = avg(x,n); s = sqrt(sum((x-m).^2/n)); end function m = avg(x,n) m = sum(x)/n; end

Function avg is a local function. Local functions are only available to other functions within the same file. Call function stat2 from the command line.

values = [12.7, 45.4, 98.9, 26.6, 53.1]; [ave,stdev] = stat2(values)

2. Adicione uma estrutura “if” para executar código se x tem dimensão maior igual a 2 3. Implementar o comando “sum(x)” utilizando a entrutura de loop “for”

Execute no Command Window:

Exercício 3 1

2

Sistemas de controle

Função de transferência: » num=[1 1]; den=[4 2 1]; G=tf(num,den);

Formato zero-polo-ganho: » Gzpk = zpk(G); [num, den]=zp2tf (z,p,k)

Sistemas com atraso no tempo: » s = tf('s'); Gd = 2*exp(-1.5*s)/(3*s +1)

Descrição no espaço de estados: » [A,B,C,D] = tf2ss(num,den); sys = ss(A,B,C,D);

Pólos e zeros de G(s): » poles=roots(den); %ou poles2=pole(G); [Z,K]=zero(G);

Sistemas de controle

Sistemas de controle • Resposta ao degrau: num=[2 1]; den=[4 2 1]; sys=tf(num,den); step(sys); %ou t=0:0.1:20; y=step(sys,t); plot(t,y) •Degrau de amplitude A: step(A*sys); Resp. impulso: impulse(sys) • Resposta a entrada arbitrária: t=0:0.1:20; u=2*exp(-0.5.*t); y=lsim(sys,u,t); plot(t,y); y=lsim(sys,u,t,x0); %resp. a entrada u e cond. inicial x0 (sys em ss) • Resposta em frequência: bode(sys); [Gm,Pm,Wcg,Wcp] = margin(sys) • Lugar das raízes: rlocus(sys) • Atraso no tempo: s=tf('s'); G = (2*exp(-0.52*s))/(10*s+1); [n d]=pade(0.52,2); %aproximação de Padé de ordem 2

• Outros: zpkdata; pzmap; rlocus, bode, margin; ltiview; c2d; d2c; printsys(num, den,‘s’)

Exercício 4 Seja

1. Expanda F(s) em frações parciais e determine a transformada inversa de Laplace. 2. Ache os pólos, zeros e o ganho de F(s). 3. Plote a resposta ao degrau de amplitude 2 de F(s). 4. Plote a resposta de F(s) à entrada r(t)=exp(-0.7*t). 5. Plote o diagrama de Bode e encontre as margens de fase e de ganho de F(s). 6. Plote o Lugar Geométrico das Raízes de F(s). 7. Forneça a representação no espaço de estados de F(s). 8. Para a representação no espaço de estados plote a resposta à condição inicial

x0=[1 2 3 4]’ 9. Para a representação no espaço de estados plote a resposta ao degrau unitário e

a x0=[1 2 3 4]’.

Conexões de modelos Comando Descrição

feedback Feedback connection of two models

connect Block diagram interconnections of dynamic systems

sumblk Summing junction for name-based interconnections

series Series connection of two models

parallel Parallel connection of two models

append Group models by appending their inputs and outputs

blkdiag Block-diagonal concatenation of models

imp2exp Convert implicit linear relationship to explicit input-output relation

inv Invert models

lft Generalized feedback interconnection of two models (Redheffer star product)

Exemplo: G = (2*exp(-0.5*s))/(s^2 + s + 1); Gpi = 0.1/s; Gma = Gpi*G; Gmf = feedback(Gma,1); step(Gmf)

Exercício 5 1. Plotar a resposta ao degrau de

𝐺 𝑠 =3𝑒−2𝑠

10𝑠 + 1

2. Seja as funções de transferência de um processo G(s), um transmissor

H(s) e um controlador de realimentação C(s)

𝐺 𝑠 =3

10𝑠+1 , H 𝑠 =

1

2𝑠+1 , C 𝑠 = 0.2( 1 +

1

5𝑠)

Forneça 1. Resposta ao degrau de amplitude 2 do sistema em malha fechada 2. Resposta a entrada u=2*sin(-0.5.*t), 0<t<20, do sistema em malha

fechada 3. Resposta do sistema em malha fechada a condição inicial x0 = 15

(sistema descrito no espaço de estados). Comandos recomendados: • tf, feedback, step, lsim, tf2ss

Exercício 6 1. Let

2. Para apresente a resposta ao degrau do sistema de malha fechada.

Diagrama de blocos (Matlab)

-> Implemente a resposta ao degrau do sistema ao lado. Ref.: Seborg, D., T. Edgar, and D. Mellichamp, Process Dynamics and Control, Wiley, 1989.

Diagrama de blocos (Matlab)

• blksys = append(sys1,sys2,...,sysN)

• sysc = connect(blksys,connections,inputs,outputs)

Ex.: i-ésima linha de ‘connections’: [7 2 -15 6]

indica que y(2)-y(15)+y(6) alimenta u(7)

Diagrama de blocos (Matlab)

sys=append(sys1,sys2,sys3,sys4,sys5,sys6,sys7,sys8,sys9);

Q=[2 1 -9; 3 2 0; 4 3 0; 5 4 0; 7 5 6; 8 7 0; 9 8 0]

inputs=[1 6]

outputs=7

sysc = connect(sys,Q,inputs,outputs);

sysr = minreal(sysc);

step(sysr)

Matrix Q: • The first line 2 1 -9 indicates

that block 2 (the PID controller) has an input which is the sum of the output from block 1 minus the output from block 9.

• The line 3 2 0 indicates that block 3 has an input from block 2.

• The other rows are also developed using the block diagram connections until all the loop elements are accounted for.

Exercício 7

1. Plote a resposta c(t) à entrada degrau unitário.

2. Obtenha a função de transferência e a resposta ao degrau de C(s)\R(s).

Simulando e.d.o’s no Matlab

Resolvendo edo’s:

• Arquivo “edo_solve.m” criado

• Se as funções estiverem em arquivos distintos parâmetros devem ser passados com argumentos de entrada ou variáveis globais

function edo_solve

[R, A, Fe] = deal(0.6,2,10); %parâmetros

h0 = 10; T = [0 200]; %cond. tnicial e tempo sim.

[T,H] = ode45(@dhdt,T,h0); %resolvendo a edo

plot(T,H); grid %plotando a edo

title('Respota cond. inicial e entrada');

xlabel('t'); ylabel('h(t)');

% Definindo a edo

function dh = dhdt(t,h)

dh = (vazao-(sqrt(h)/R))/A;

% Sinal de entrada que varia com o tempo

function v = vazao

if t < 120

v = Fe;

else

v = 0.5*Fe;

end

end

end

end

Arquivo “edo_solve.m”:

Simulando e.d.o’s no Matlab Considere a edo não-linear com entrada u:

Aplicar entrada u(t)=0.1, t < 0, e u(t)=2+0.5*sin(20*t) para t >= 10:

x0 = [0 0.25];

TSPAN = [0 40];

[t,x] = ode23('dxdt',TSPAN,x0);

% Plotando x1 e x2

plot(t,x(:,1),'r',t,x(:,2),'b');

% Plotando sinal de entrada

ti = TSPAN(1):0.01:TSPAN(2);

for i = 1:length(ti)

if ti(i) > 10

u(i)=2+0.5*sin(20*ti(i));

else u(i)=0.1; end

end

figure; plot(ti,u);

function dx = dxdt(t,x)

if t > 10, u=2+0.5*sin(20*t);

else u=0.1; end

dx1 = -x(1)*(1-x(2)^2)-x(2);

dx2 = x(1) + u;

dx = [dx1 dx2]';

Criar arquivo “dxdt.m”:

Criar e executar arquivo “edononlinear.m”:

Simulando e.d.o’s no Matlab

Ex: options = odeset('RelTol',1e-4,'AbsTol',1e-9);

[T,H] = ode45(@dhdt,[0 Tfinal],h0,options);

[T,X] = ode23('dxdt',[0 40],[1 2],options);

Linearização

• Exemplo:

% Sistema não-linear: dot{x} = f(x,u)

syms x1 x2 u1 u2

dx1 = -x1*(1-x2^2)-x2 + 1/u1;

dx2 = sqrt(x1) + u1*u2;

dx = [dx1; dx2];

% Matriz Jacobiana

Aj = jacobian(dx,[x1 x2]);

Bj = jacobian(dx,[u1 u2]);

% Linearização em torno de

% (x1bar,x2bar,u1bar,u2bar)= (1.2, 3.4, 5, 7)

A = eval(subs(Aj,{x1,x2,u1,u2},{1.2,3.4,5,7}));

B = eval(subs(Bj,{x1,x2,u1,u2},{1.2,3.4,5,7}));

% Sistema linear: dxtil = A*xtil + B*util

% xtil = x – xbar; util = u - ubar

sys_ss = ss(A,B,eye(2),zeros(2))

G = tf(sys_ss)

Simulink

Simulando e.d.o’s no Simulink

Biblioteca de Funções

Funções usuais

Simulando...

Exporta variáveis para o Workspace

do Matlab

Scope:

Simulando e.d.o’s no Simulink E.D.O

Diagrama Simulink:

Condições iniciais definidas no bloco integrador.

Simulando e.d.o’s no Simulink

Simulando e.d.o’s no Simulink a partir do Matlab

Processo em malha fechada

• Degrau aplicado somente após processo estar em regime permanente (idealmente no ponto de operação em t = 0)

• Condições iniciais (c.i) nos blocos de integração (c.i = ponto de operação) • Saída de uma função de transferência é nula em t=0 (variável de desvio)

• Entrada da função de transferência: variável absoluta variável de desvio • Saída da função de transferência: variável de desvio variável absoluta

• Saída do controlador • Soma-se o valor em regime permanente do sinal de controle • Em % na faixa entre 0 e 100% (bloco saturador para garantir faixa)

• Adaptação do sinal de referência em função do ganho do transmissor • Variáveis de entrada e saída dos blocos do Simulink devem ter siginificado e

unidade de engenharia conhecido

Exercício 6

Simulação no Simulink Linearização: comando jacobian Exemplo:

>> syms x y z; >> A = jacobian([x*y*z; y; x+z],[x y z]); >> Ax = eval(subs(A,{x,y,z},{1,2,3}));

Bibliografia • Mathworks, http://www.mathworks.com/help/ • Kermit Sigmon, Matlab Primer, 3ª ed,

http://www.math.toronto.edu/mpugh/primer.pdf • Rao V. Dukkipati, Matlab – An Introduction with Applications, Ed.

New Age, 2010. • Mathworks, Matlab Primer,

https://www.mathworks.com/help/pdf_doc/matlab/getstart.pdf • C. A. Smith e A. Corripio, Príncipios e Prática do Controle

Automático de Processo, 3ª. Edição, Ed. LTC, 2012. • Robinson, T.; Kambouchev, N., 16.06/16.07 Matlab/Simulink

Tutorial, Massachusetts Institute of Technology, http://dspace.mit.edu/bitstream/handle/1721.1/60691/16-07-fall-2004/contents/study-materials/matlabtut.pdf

• Apostila de Introdução ao MATLAB, Universidade Federal Fluminense, Centro Tecnológico, Escola de Engenharia, http://www.telecom.uff.br/pet/petws/downloads/apostilas/MATLAB.pdf

• Stephen J. Chapman, Essentials of MATLAB Programming, 2nd ed. • Michael G. Kay, Basic Concepts in Matlab, 2017