%% Cancelaciones en un controlador IMC (discreto), relación con controlabilidad/observabilidad de un sistema en bucle cerrado
% © 2021, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Presentación en vídeo: http://personales.upv.es/asala/YT/V/cancobsctr.html
% Este código funcionó correctamente con Matlab R2021b
% Objetivo: Las cancelaciones polo-cero entre distintos elementos (en este ejemplo, entre plantas y controladores) hacen cosas no observables o no controlables desde determinadas entradas y salidas... pero desde otras quizás no contempladas (o condiciones iniciales no nulas), sí pueden verse sus efectos.
% Table of Contenidos
% Modelo
% Diseño (no recomendable) de control IMC discreto
% Prestaciones en bucle cerrado
% Controlabilidad y observabilidad del bucle cerrado
%     Forma canónica modal
% Conclusión
%% Modelo

Ac=[0 1 ;-3 -0.2]; Bc=[0;3]; C=[1 0];
sysc=ss(Ac,Bc,C,0);sysc.InputName='ureal';sysc.OutputName='y';
tf(sysc) %vemos la TF, pero operaremos en repr. interna.

% Vamos a discretizarlo:

Ts=0.2; 
sysd=c2d(sysc,Ts,'zoh');
zpk(sysd)

%
%% Diseño (no recomendable) de control IMC discreto

s=tf('s');
ModeloRefContinuo=6/(s^2+3.6*s+6);
ModeloRef=c2d(ModeloRefContinuo,Ts,'zoh');
zpk(ModeloRef)
step(sysd,ModeloRef,40), grid on, legend("Bucle abierto", "Modelo Ref.")
Q=minreal(ModeloRef/sysd); %parámetro de Youla del control por modelo interno

% El regulador IMC será:

K_IMC=minreal(feedback(Q,-sysd)); size(K_IMC)

%% Prestaciones en bucle cerrado
%
% Cerramos bucle con "connect", ante perturbación de entrada.

K_IMC.InputName='e';
K_IMC.OutputName='u';
s1=sumblk('e=r-y');s2=sumblk('ureal=u+du');
BucleCerrado=connect(s1,s2,K_IMC,sysd,{'r','du'},{'y','u'});
step(BucleCerrado,40), grid on
bodemag(BucleCerrado), grid on

%
% No hay buena respuesta por la cancelación planta/controlador:

zpk(sysd)
zpk(K_IMC)

%
%% Controlabilidad y observabilidad del bucle cerrado
% ¿Se podía detectar estas cosas raras que ocurren en el step/bodemag basándose en controlabilidad/observabilidad?

size(BucleCerrado)
hsvd(BucleCerrado)' %todo es "controlable" y "observable"
svd(ctrb(BucleCerrado.A,BucleCerrado.B))' %todo es "controlable"
rank(ctrb(BucleCerrado.A,BucleCerrado.B))' %todo es "controlable"
svd(obsv(BucleCerrado.A,BucleCerrado.C))' %todo es "observable"
rank(obsv(BucleCerrado.A,BucleCerrado.C))' %todo es "observable"

% El sistema sea "controlable" y "observable" (esto es, todo lo que pasa por dentro lo mueve "r o d_u", y se ve en "y o u")... o sea, que lo que hemos visto en step/bodemag es "todo lo que hay que ver": si algo raro pasa internamente en el sistema, en alguna de las 4 gráficas se ha de ver: no necesitamos un diagrama de bloques más complicado con otras entradas o salidas para estudiar este bucle.
% Pero, podemos sacar más conclusiones usando esto de controlabilidad y observabilidad, para entender en más detalle lo que está pasando.

svd(ctrb(BucleCerrado.A,BucleCerrado.B(:,1)))' %2 estados no son "controlables" desde "r"
svd(ctrb(BucleCerrado.A,BucleCerrado.B(:,2)))' %1 estado no es "controlable" desde "du"
svd(obsv(BucleCerrado.A,BucleCerrado.C(1,:)))' %1 estado no es "observable" desde "y"
svd(obsv(BucleCerrado.A,BucleCerrado.C(2,:)))' %2 estados no son "observables" desde "u"

% Esto indica que si sólo nos fijaramos en una de las entradas o en una de las salidas, nos dejaríamos algo "en el tintero" sobre la dinámica interna de este bucle.
% Para ver qué estados (o combinaciones de ellos) no son controlables u observables desde qué entradas o salidas, pasemos a forma modal.
%% Forma canónica modal

eBC=eig(BucleCerrado)'
abs(eBC)

% Pasando a forma canónica modal (diagonalizando; bueno, los pares de polos complejos de 2 en 2 con [re im;-im re]):

[Canon_sys, Ti]=canon(BucleCerrado,'modal');

% A priori, los estados pueden ser mezcla de estados de planta y controlador... no supone inconveniente, porque las entradas y salidas no cambian:

Ti

% Esta es la forma canónica final:

Canon_sys

%
% Las condiciones iniciales en cada uno de los nuevos estados dan lugar a los siguientes transitorios:

initial(Canon_sys,[1;0;0;0;0],40) %condiciones iniciales coincidentes con autovectores ( identidad en nuevas coordenadas)

% El polo se ve en "u" pero no en "y".
%

initial(Canon_sys,[0;1;0;0;0],40) %condiciones iniciales coincidentes con autovectores ( identidad en nuevas coordenadas)

% El polo se ve en "y" pero no en "u".
%

initial(Canon_sys,[0;0;0;1;0],40) %condiciones iniciales coincidentes con autovectores ( identidad en nuevas coordenadas)

% Esta dinámica se ve tanto en "y" como en "u".
%
%% Conclusión
%
% \bulletHay dos modos (polos) que tardan 40 segundos... y que no se hace control para amortiguarlos... Los componentes de condiciones iniciales en esas direcciones darán una dinámica mala, también si la excitan los ruidos de proceso (pert. entrada): al final de los transitorios sólo quedarán esos modos... la referencia no los excita, pero las condiciones inciales fuera de equilibrio, posibles perturbaciones y errores de modelado que siempre habrá en la práctica, es posible que sí.
% \bulletHay otro modo no observable en la salida, pero que da lugar a una acción de control de alta frecuencia. Tampoco es deseable.
% Esto justifica las reglas usuales "de libro" sobre IMC o cualquier otro regulador basado en cancelación (p.ej. mixed sensitivity): (a) No uses IMC en procesos para los que desees cambiar "mucho" la dinámica de bucle abierto porque cambiará ante referencia, pero no ante "initial" o "perturb. entrada"... no canceles polos cerca del límite de inestabilidad; (b) No canceles ceros cerca del límite de inestabilidad, porque te destrozará la respuesta en acción de control.