Aula 7 - - ITA, Divisão de Engenharia Eletrônicakawakami/ee254/EE254_2oSem_2019_Aula7...EE-254 (Controle Preditivo) Aula 7 10 Setembro 2019 22 / 59 Convexidade Para mais detalhes,

Aula 7

10 Setembro 2019

EE-254 (Controle Preditivo) Aula 7 10 Setembro 2019 1 / 59

Resumo da aula passada

Estimacao de perturbacoes empregando observador de estados

Determinacao de valores de equilıbrio para o estado e o controle


Topicos da aula de hoje

Tratamento de restricoes - caso SISO


Restricoes a serem consideradas

Tres tipos basicos de restricoes serao considerados:

Incrementos no controle ∆u

Excursao do controle u

Excursao da saıda y

Como se vera, todas essas restricoes podem ser expressas em termos derestricoes sobre ∆u.


Restricoes sobre os incrementos no controle ∆u

∆umin ≤ ∆u(k + i − 1|k) ≤ ∆umax , i = 1, 2, . . . ,M

Vale ressaltar que ∆umin e ∆umax sao limitantes para o incremento nocontrole por perıodo de amostragem.

Ex: Suponha que o controle u corresponda a uma deflexao medida emgraus e que a taxa de variacao de u esteja limitada a ± 10◦/ s.

Se o perıodo de amostragem for T = 10 ms, quais serao os valores de∆umin e ∆umax ?

Resp: ∆umin = −0,1◦ e ∆umax = +0,1◦.


Notacao alternativa:

1M︷︸︸︷11...1

∆umin ≤

∆u︷︸︸︷∆u(k|k)

∆u(k + 1|k)...

∆u(k + M − 1|k)

≤1M︷︸︸︷11...1

∆umax

1M∆umin ≤ ∆u ≤ 1M∆umax

Vale ressaltar que o sımbolo ≤ e aqui empregado para denotar umadesigualdade elemento a elemento.


1M∆umin ≤ ∆u ≤ 1M∆umax{∆u ≤ 1M∆umax

−∆u ≤ −1M∆umin

2M×M︷︸︸︷[IM−IM

]∆u ≤

2M×1︷︸︸︷[1M∆umax

−1M∆umin

]


Restricoes sobre a excursao do controle u

umin ≤ u(k + i − 1|k) ≤ umax , i = 1, 2, . . . ,M

Obs: Se o modelo tiver sido linearizado em torno de um valor de equilıbrioup para o controle (vide Aula 3), os limitantes umax e umin corresponderaoa diferencas com respeito a up.

Ex: Suponha que o controle seja gerado por uma fonte de tensao quesatura em up,min = −5 V e up,max = +5 V.

Se a operacao se der em torno de um valor de equilıbrio up = +3 V, quaisserao os valores de umin e umax a serem empregados nas restricoes ?

Resp: umin = up,min − up = −8 V e umax = up,max − up = +2 V.



1Mumin ≤ u ≤ 1Mumax

Para expressar estas restricoes em termos de ∆u, deve-se obter umarelacao entre u e ∆u.


u(k|k) = u(k − 1) + ∆u(k|k)

u(k + 1|k) = u(k|k) + ∆u(k + 1|k)

= u(k − 1) + ∆u(k|k) + ∆u(k + 1|k)...

u(k + M − 1|k) = u(k − 1) + ∆u(k|k) + · · ·+ ∆u(k + M − 1|k)


u(k|k)

u(k + 1|k)...

u(k + M − 1|k)

=

u(k − 1)u(k − 1)

...u(k − 1)

+

+

∆u(k|k)∆u(k|k) + ∆u(k + 1|k)...∆u(k|k) + ∆u(k + 1|k) + · · ·+ ∆u(k + M − 1|k)


u(M×1)︷︸︸︷u(k|k)

u(k + 1|k)...

u(k + M − 1|k)

=

1Mu(k−1)︷︸︸︷u(k − 1)u(k − 1)

...u(k − 1)

+

1 0 · · · 01 1 · · · 0...

.... . .

...1 1 · · · 1

︸︷︷︸

TM(M×M)

∆u(k|k)

∆u(k + 1|k)...

∆u(k + M − 1|k)

︸︷︷︸

∆u(M×1)

u = 1Mu(k − 1) + TM∆u


1Mumin ≤ u ≤ 1Mumax

1Mumin ≤ 1Mu(k − 1) + TM∆u ≤ 1Mumax

1M [umin − u(k − 1)] ≤ TM∆u ≤ 1M [umax − u(k − 1)]

{TM∆u ≤ 1M [umax − u(k − 1)]

−TM∆u ≤ 1M [u(k − 1)− umin]

2M×M︷︸︸︷[TM

−TM

]∆u ≤

2M×1︷︸︸︷[1M [umax − u(k − 1)]1M [u(k − 1)− umin]

]


Equıvoco a ser evitado

Uma restricao da forma

TM∆u ≤ 1M [umax − u(k − 1)]

nao e equivalente a

∆u ≤ T−1M 1M [umax − u(k − 1)]


Exemplo: M = 2, u(k − 1) = 0, umax = 1

TM =

[1 01 1

], T−1

M =

[1 0−1 1

]Forma correta de se escrever a restricao:

TM∆u ≤ 1M [umax − u(k − 1)][1 01 1

] [∆u(k|k)

∆u(k + 1|k)

]≤[

11

]{

∆u(k|k) ≤ 1

∆u(k|k) + ∆u(k + 1|k) ≤ 1


Exemplo: M = 2, u(k − 1) = 0, umax = 1

Forma correta de se escrever a restricao:{∆u(k|k) ≤ 1

∆u(k|k) + ∆u(k + 1|k) ≤ 1

Forma incorreta de se escrever a restricao:

∆u ≤ T−1M 1M [umax − u(k − 1)][

∆u(k|k)∆u(k + 1|k)

]≤[

1 0−1 1

] [11

]{

∆u(k|k) ≤ 1

∆u(k + 1|k) ≤ 0


Forma correta de se escrever a restricao:{∆u(k|k) ≤ 1∆u(k|k) + ∆u(k + 1|k) ≤ 1


Forma incorreta de se escrever a restricao:{∆u(k|k) ≤ 1∆u(k + 1|k) ≤ 0


Restricoes sobre a excursao da saıda y

ymin ≤ y(k + i |k) ≤ ymax , i = 1, 2, . . . ,N

Obs: Se o modelo tiver sido linearizado em torno de um valor de equilıbrioyp para a saıda (vide Aula 3), os limitantes ymax e ymin corresponderao adiferencas com respeito a yp, isto e

ymin = yp,min − yp

ymax = yp,max − yp



1Nymin ≤ y ≤ 1Nymax

Lembrando que y = G∆u + f, a restricao pode ser reescrita como

1Nymin ≤ G∆u + f ≤ 1Nymax

1Nymin − f ≤ G∆u ≤ 1Nymax − f

{G∆u ≤ 1Nymax − f

−G∆u ≤ f − 1Nymin

2N×M︷︸︸︷[G−G

]∆u ≤

2N×1︷︸︸︷[1Nymax − ff − 1Nymin

]EE-254 (Controle Preditivo) Aula 7 10 Setembro 2019 20 / 59

Resumo das restricoes

S[(4M+2N)×M]︷︸︸︷

IM−IMTM

−TM

G−G

∆u ≤

b[(4M+2N)×1]︷︸︸︷

1M∆umax

−1M∆umin

1M [umax − u(k − 1)]1M [u(k − 1)− umin]

1Nymax − ff − 1Nymin

S∆u ≤ b


Formulacao do problema de otimizacao com restricoes

O problema de otimizacao na presenca das restricoes consideradas podeser formulado como

min∆u∈RM

J(∆u) =1

2∆uTH∆u + cT∆u + cte

sujeito aS∆u ≤ b

Funcao de custo quadratica com H > 0

Restricoes lineares

Problema de Programacao Quadratica

Pode-se mostrar que este e um problema de otimizacao convexo.


Convexidade

Para mais detalhes, vide Capıtulos 2, 3 e 4 da seguinte referencia:

Boyd, S.; Vandenberghe, L. Convex Optimization. Cambridge UniversityPress, 2004.


Conjuntos convexos

Um conjunto C e dito ser convexo se, para quaisquer x , y ∈ C e qualquerescalar λ tal que 0 ≤ λ ≤ 1, tem-se

λx + (1− λ)y ∈ C


Conjunto definido por interseccao de semiespacos

Pode-se mostrar que um conjunto C nao vazio definido como

C ={x ∈ RM

∣∣∣ Sx ≤ b}

e convexo.

Com efeito, dados x , y ∈ C e λ ∈ [0, 1], tem-se

S [λx + (1− λ)y ] = λSx + (1− λ)Sy ≤ λb + (1− λ)b = b

Portanto λx + (1− λ)y ∈ C.


Funcoes convexas

Uma funcao f : C → R e dita ser convexa se C for um conjunto convexo e

f(λx + (1− λ)y

)≤ λf (x) + (1− λ)f (y)

para quaisquer x , y ∈ C e λ ∈ [0, 1].

A funcao f e dita ser estritamente convexa se

f(λx + (1− λ)y

)< λf (x) + (1− λ)f (y)

sempre que x 6= y e 0 < λ < 1.


Funcao quadratica

Pode-se mostrar que uma funcao f : C → R da forma

f (x) = xTQx + cT x + d , Q = QT > 0

e estritamente convexa. Com efeito, sejam x , y ∈ C e λ ∈ R, com x 6= y e0 < λ < 1. Tem-se, entao:

f(λx + (1− λ)y

)=

(λx + (1− λ)y

)TQ(λx + (1− λ)y

)+ cT

(λx + (1− λ)y

)+ d

= λ2xTQx + λ(1− λ)xTQy + λ(1− λ)yTQx + (1− λ)2yTQy

+λcT x + (1− λ)cT y + d


Por outro lado:

λf (x) + (1− λ)f (y) =

λ(xTQx + cT x + d

)+ (1− λ)

(yTQy + cT y + d

)=

λxTQx + (1− λ)yTQy + λcT x + (1− λ)cT y + d


Logo:

f(λx + (1− λ)y

)−[λf (x) + (1− λ)f (y)

]=

λ2xTQx + λ(1− λ)xTQy + λ(1− λ)yTQx + (1− λ)2yTQy

(((((((((((((+ λcT x + (1− λ)cT y + d

−λxTQx − (1− λ)yTQy(((((((((((((− λcT x − (1− λ)cT y − d

= λ(λ− 1)xTQx − λ(λ− 1)xTQy − λ(λ− 1)yTQx + (1− λ)(−λ)yTQy

= λ(λ− 1)[(x − y)TQ(x − y)

] ∗< 0

*Q > 0, x 6= y , 0 < λ < 1


Mınimo local de uma funcao convexa

Seja C um conjunto convexo e f : C → R uma funcao convexa. Suponhaque x ′ ∈ C seja um ponto de mınimo local de f , isto e:

f (x ′) = minx∈C, ||x−x ′||≤r

f (x)

para algum r > 0. Entao, pode-se mostrar que x ′ tambem e um ponto demınimo global.


f (x ′) = minx∈C, ||x−x ′||≤r

f (x) (1)

Por absurdo, suponha que exista x ′′ ∈ C tal que

f (x ′′) < f (x ′)

Nesse caso, em vista de (1), tem-se

||x ′′ − x ′|| > r > 0 (2)

Seja agora y = (1− λ)x ′ + λx ′′ com

λ =r

2||x ′′ − x ′||(3)

De (2) e (3), segue que 0 < λ < 1.


y = (1− λ)x ′ + λx ′′

λ =r

2||x ′′ − x ′||

Como C e convexo e 0 < λ < 1, tem-se que y ∈ C. Adicionalmente:

||y − x ′|| = || − λx ′ + λx ′′|| = λ||x ′′ − x ′|| =r

2< r


f (x ′) = minx∈C, ||x−x ′||≤r

f (x) (1)

||y − x ′|| < r (4)

Como y = (1− λ)x ′ + λx ′′, com 0 < λ < 1, e a funcao f e convexa,tem-se que

f (y) ≤ (1− λ)f (x ′) + λ f (x ′′)︸︷︷︸<f (x ′)

< f (x ′)

o que contradiz (1) e (4).


Mınimo local de uma funcao estritamente convexa

Seja C um conjunto convexo e f : C → R uma funcao estritamenteconvexa. Suponha que x ′ ∈ C seja um ponto de mınimo local de f , isto e:

f (x ′) = minx∈C, ||x−x ′||≤r

f (x)

para algum r > 0. Entao, se x ∈ C e f (x) = f (x ′), pode-se mostrar quex = x ′.


Se x ∈ C e f (x) = f (x ′), pode-se mostrar que x = x ′.

Por absurdo, suponha que x 6= x ′. Seja entao

y = λx + (1− λ)x ′

com 0 < λ < 1. Da convexidade de C, tem-se que y ∈ C. Sabendo que f eestritamente convexa, conclui-se que

f (y) < λ f (x)︸︷︷︸f (x ′)

+(1− λ)f (x ′) = f (x ′)

o que contradiz a hipotese de que x ′ e mınimo local (e, portanto, global)da funcao f .


Programacao Quadratica (Ex: M = 2)


Programacao Quadratica no Matlab

Funcao QUADPROG (Optimization Toolbox):

x = quadprog(H,f,A,b)

minimiza 0.5*x’*H*x + f’*x

sujeito a A*x <= b

Em nosso caso:

xqp = ∆uHqp = H = 2(GTG + ρIM)fqp = c = 2GT (f − r)Aqp = Sbqp = b

podendo ser omitido o fator 2 em H e f.


Argumentos adicionais da funcao QUADPROG

X = QUADPROG(H,f,A,b,Aeq,beq,LB,UB,X0,options)

X = QUADPROG(H,f,A,b,Aeq,beq) solves the problem above

while additionally satisfying the equality constraints

Aeq*x = beq.

X = QUADPROG(H,f,A,b,Aeq,beq,LB,UB) defines a set of

lower and upper bounds on the design variables, X, so

that the solution is in the range LB <= X <= UB. Use

empty matrices for LB and UB if no bounds exist. Set

LB(i) = -Inf if X(i) is unbounded below; set UB(i) =

Inf if X(i) is unbounded above.

X0: starting point.

options: opcoes de otimizacao, incluindo algoritmo a ser usado ecriterios de parada.


Opcoes de algoritmos:

trust-region-reflective (formerly LargeScale = ’on’): default

active-set (formerly LargeScale = ’off’)

interior-point-convex


Programacao Quadratica no Matlab: Exemplo

minx1,x2

J =1

2

[(x1 − 5)2 + (x2 − 5)2

]s.a.

0 ≤ x1 ≤ 2

0 ≤ x2 ≤ 2



Custo:

J =1

2

[(x1 − 5)2 + (x2 − 5)2

]=

1

2(x2

1 − 10x1 + 25 + x22 − 10x2 + 25)

=1

2(x2

1 + x22 )− 5x1 − 5x2 + 25

1

2[x1 x2]

[1 00 1

]︸︷︷︸

Hqp

[x1

x2

]+ [−5 − 5]︸︷︷︸

fqpT

[x1

x2

]+ 25


Restricoes:

0 ≤ x1 ≤ 2

0 ≤ x2 ≤ 2

x1 ≤ 2

x2 ≤ 2

−x1 ≤ 0

−x2 ≤ 0

1 00 1−1 00 −1

︸︷︷︸

Aqp

[x1

x2

]≤

2200

︸︷︷︸

bqp


>> H = eye(2);

>> f = [-5;-5];

>> b = [2*ones(2,1);0*ones(2,1)];

>> A = [eye(2);-eye(2)];

>> x = quadprog(H,f,A,b);

Resultado: x = [2;2].


Para especificar o algoritmo a ser empregado:

>> options = optimset(‘Algorithm’,‘active-set’);

>> quadprog(H,f,A,b,[],[],[],[],[],options);


Programacao Quadratica: Alguns pacotes computacionais

CVX (Matlab software for disciplined convex programming):cvxr.com/cvx

CVXGEN (Geracao de codigo em C): cvxgen.com

IPOPT (Codigo em C++ empregando metodo de pontos interiores):www.coin-or.org/Ipopt

QL (Codigo original em Fortran, traduzido posteriormente para C):www.ai7.uni-bayreuth.de/software.htm

CPLEX (IBM):www-01.ibm.com/software/integration/optimization/cplex-optimizer/

GUROBI: www.gurobi.com


MPC no espaco de estados com restricoes

Informacao requerida sobre a planta:

Matrizes A,B,C do modelo no espaco de estados

Limitantes sobre os incrementos no controle: ∆umin,∆umax

Limitantes sobre a excursao do controle: umin, umax

Limitantes sobre a excursao da saıda: ymin, ymax

Parametros de projeto:

Peso do controle ρ

Horizonte de predicao N

Horizonte de controle M


Inicializacao:

Fazer

A =

[A B

0Tn 1

], B =

[B1

], C = [C 0]

G =

C B 0 · · · 0

C AB C B · · · 0...

.... . .

...

C AN−1B C AN−2B · · · C AN−M B

, Φ =

C A

C A2

...

C AN


Fazer Hqp = 2(GTG + ρI ) , Aqp =

IM−IMTM

−TM

G−G

Fazer k = 0, u(−1) = 0


Rotina principal:

1 Ler x(k) (estado da planta) e yref (valor de referencia para a saıda)

2 Fazer r = [yref ]N3 Fazer

ξ(k) =

[x(k)

u(k − 1)

]4 Calcular f = Φ ξ(k) e fqp = 2GT (f − r)

5 Fazer bqp =

1M∆umax

−1M∆umin

1M [umax − u(k − 1)]1M [u(k − 1)− umin]

1Nymax − ff − 1Nymin


6 Resolver o problema de otimizacao

∆u∗ = arg min∆u∈RM

1

2∆uTHqp∆u + f Tqp∆u s.a. Aqp∆u ≤ bqp

7 Obter o incremento no controle:

∆u(k) = ∆u∗(k|k) = [ 1 0 · · · 0 ] ∆u∗

8 Atualizar o controle aplicado a planta: u(k) = u(k − 1) + ∆u(k)

9 Fazer k = k + 1

10 Aguardar o proximo instante de amostragem e retornar ao passo 1.


Implementacao em Matlab

matrizes_ss_du_restricoes.m: Monta as matrizes Φ,G ,Hqp,Aqp.

mpc_ss_du_restricoes.m: S-function que implementa o controlador


Exemplo: Sistema de levitacao magnetica

Arquivos Matlab:

parametros_maglev_ss.m: Define os parametros do levitadormagnetico

levitador_mpc_ss_du_restricoes.mdl: Diagrama de simulacao

levitador_mpc_ss_du_restricoes_pert.mdl: Diagrama desimulacao com estimacao de perturbacoes


Importancia do tratamento de restricoes no controle

A importancia de tratar restricoes na saıda y esta usualmente relacionadaa imposicoes de seguranca ou requisitos de qualidade para o produto final.

Qual a necessidade de tratar restricoes sobre a entrada (em termos deexcursao e/ou incrementos) ?

Por que nao calcular a sequencia de controle (ou incrementos decontrole) otima e satura-la nos limitantes inferior e superior, senecessario ?




Restricoes no controle: Acomodacao de falhas

Em particular, a possibilidade de tratar restricoes no controle pode serutil para acomodacao de falhas de atuadores.

Nesse contexto, um sistema de monitoramento de falhas informaria aocontrolador preditivo os novos limitantes para a excursao e/ouincrementos do sinal de controle.

Exemplo: levitador_mpc_ss_du_restricoes_falha.mdl


Resumo da aula de hoje

Tratamento de restricoes - Caso SISO:

Incrementos no controle ∆u

Excursao do controle u

Excursao da saıda y

Programacao quadratica

Convexidade

Implementacao em Matlab


Proxima aula: Prova

Duracao: 150 minutos

Consulta permitida a livros, anotacoes pessoais e material distribuıdoao longo do curso.


Documents

Aula 7 - - ITA, Divisão de Engenharia Eletrônicakawakami/ee254/EE254_2oSem_2019_Aula7...EE-254 (Controle Preditivo) Aula 7 10 Setembro 2019 22 / 59 Convexidade Para mais detalhes,