%% Modelado para control robusto: alternativas con distinto numero de incidencias de un párametro incierto
% © 2020, 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/urealmult.html
% Este código funcionó correctamente con Matlab R2020a
% Objetivo: diferentes formas de escribir las ecuaciones de un modelo dinámico con parámetros inciertos pueden llevar a realizaciones no mínimas (que solventaríamos con minreal) pero también a diferente número de veces (incidencias) que aparece un determinado parámetro incierto. El comando simplify de la Robust Control Toolbox intenta reducir dicho número de incidencias, pero todo resulta más fácil si las ecuaciones de la física son agrupadas de modo que el parámetro incierto aparece el menor número de veces posible.
% *Nota: en versiones posteriores de Matlab, es posible que sea más "inteligente" con el tema multiincidencia y salgan menos "occurrences" que en este mlx/PDF guardado.
% Tabla de Contenidos
% Modelo (simbólico) de un péndulo linealizado
%     Sustitución de valores de parámetros (conocidos o inciertos)
% Alternativa: forma "manual" a partir de diagrama de bloques linealizado
% Comparación de prestaciones/tiempo de cómputo con la Robust Control Toolbox
% Conclusiones
%
%% Modelo (simbólico) de un péndulo linealizado

syms dx1dt dx2dt x1 x2 M L F g b
modelo=[dx1dt==x2;M*L^2*dx2dt==-sin(x1)*L*M*g-b*x2+F]
EcEstado=solve(modelo,{dx1dt,dx2dt});
dxdt=[EcEstado.dx1dt; EcEstado.dx2dt]

% Para linealizar la ecuación de estado hacemos derivadas parciales:

Asym=jacobian(dxdt,[x1 x2])
Bsym=jacobian(dxdt,F)

% Sustituimos ahora por valores en punto de equilibrio:

x1eq=0;
Asym_eq=subs(Asym,x1,x1eq)

% En B no hace falta, porque no sale x1.
%% Sustitución de valores de parámetros (conocidos o inciertos)
% Hacemos los parámetros constantes o inciertos:

g=9.81; 
b=ureal("b",0.7,"PlusMinus",[-.1,.1]);b.AutoSimplify='full';
L=ureal("L",1.25,"Percent",2);L.AutoSimplify='full';
M=ureal("M",0.5,"Percent",5);M.AutoSimplify='full';

% Nota importante: estamos incorporando la incertidumbre paramétrica, pero no la debida al error de linealización.
% Los sustituimos en las expresiones simbólicas para obtener objetos de la robust control toolbox:

A=evalin('base',char(Asym_eq))
A.NominalValue
B=evalin('base',char(Bsym))
B.NominalValue
sys=ss(A,B,[1 0],0) %escogemos posición angular como salida
step(sys,8), grid on

% Intentamos reducir (sin éxito) el número de incidencias de alguno de los parámetros:

simplify(sys,'full')

%
%
%% Alternativa: forma "manual" a partir de diagrama de bloques linealizado
% Podemos considerar este diagrama de bloques, con 2 incidencias de "M", y 3 de "L", y una de "b":
%
% Pero podemos transformarlo en este otro con 1 incidencia de "M", dos de "L" y 1 de "b":
%
%

s=tf('s');
integrador1=1/s; integrador1.InputName={'x2'};integrador1.OutputName={'x1'};
integrador2=1/s; integrador2.InputName={'r'};integrador2.OutputName={'x2'};
unoentreele1=tf(1/L);unoentreele1.InputName={'q'};unoentreele1.OutputName={'r'};
sum1=sumblk('q=n-m');
unoentreele2=tf(1/L);unoentreele2.InputName={'t'};unoentreele2.OutputName={'n'};
gblock=tf(g);gblock.InputName={'x1'};gblock.OutputName={'m'};
unoentreeme=tf(1/M);unoentreeme.InputName={'a'};unoentreeme.OutputName={'t'};
sum2=sumblk('a=F-w');
bblock=tf(b);bblock.InputName={'x2'};bblock.OutputName={'w'};
sys2=connect(integrador1,integrador2,unoentreele1,unoentreele2,sum1,gblock,unoentreeme,sum2,bblock,'F','x1')
norm(sys2.NominalValue-sys.NominalValue,'inf')

%% Comparación de prestaciones/tiempo de cómputo con la Robust Control Toolbox

opt = robOptions('Display','on','Sensitivity','on');
tic,[rs,wcus]=robstab(sys,opt),toc
tic,[rs,wcus]=robstab(sys2,opt),toc

%% Conclusiones
% Repeticiones innecesarias de los parámetros inciertos causan mayor tiempo de cómputo y resultados más conservativos cuando la Robust Control Toolbox se usa para analizar estabilidad robusta o diseñar controladores (musyn). Es, pues, muy recomendable evitar dichas repeticiones innecesarias en la fase de modelado.
% Recalcamos de nuevo que la incertidumbre no lineal no ha sido modelada: en la documentación de Matlab existe el bloque udyn de incertidumbre no-lineal, variante en el tiempo, pero las rutinas de incertidumbre estructurada como robstab no lo manejan.