%% Reducción de orden: influencia de modos longitudinales en comportamiento de muelles
% © 2019 Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Objetivo: comprender cómo simplificar un modelo mecánico de alto orden (elementos finitos) de modo que se pierda poca exactitud para su uso posterior.
% Presentación en vídeo en: http://personales.upv.es/asala/YT/V/feredm.html
% http://personales.upv.es/asala/YT/V/feredm2.html
%
% Table of Contents
% 1.-Planteamiento del problema y modelado
% 2.-Estudio numérico de un ejemplo concreto
%         *Utilización (formalmente incorrecta) de minreal
%     2.1-Respuesta temporal/frecuencia del modelo de muelle (sin carga adicional)
%     2.2-Reducción basada en realización equilibrada del modelo de orden 80 (40 elementos)
%     2.3-Reducción basada en la separación de frecuencias (modal)
% 3.-Utilización del modelo reducido del muelle como subsistema
% Apéndice: funciones auxiliares
%% 1.-Planteamiento del problema y modelado
% Un muelle no es sólo F=k(x-L), tiene su propia masa y estructura interna.
% En concreto, puede presentar modos de vibración longitudinales:
%
% ¿Influye esto en las aplicaciones?
% Modelemos el muelle, como N segmentos de muelle:
%
% Este tipo de metodologías son la idea que originó el enfoque de "elementos finitos" en mecánica, fluidos, calor, etc. Cuando los elementos son infinitesimales, se llega a una ecuación en derivadas parciales (si no hay rozamiento, a la ecuación de onda unidimensional,https://en.wikipedia.org/wiki/Wave_equation )
% \frac{\partial^2 u}{\partial t^2}=\frac{KL^2}{M}\cdot \frac{\partial^2 u}{\partial x^2}
% donde u(x,t) es el desplazamiento incremental sobre su posición natural del elemento que en el muelle sin carga estaría en la posición x.
% Las ecuaciones (linealizadas, eliminando longitud natural) de la cadena de N masas-muelle-amortiguador está modelada en la función auxiliar modelobarra al final de este fichero.
% El orden del sistema es 2N: las variables de estado son N posiciones y N velocidades.
%% 2.-Estudio numérico de un ejemplo concreto

Mtotal=.1;
Ktotal=250;  %constante elástica de la barra.
friccion_tot=0.4; %fricción

% Modelado como 40 elementos:

barra40=modelobarra(40,Ktotal,Mtotal,friccion_tot);
size(barra40)

%% *Utilización (formalmente incorrecta) de minreal
% Importante: minreal es para corregir modelos mal planteados o mal manipulados. Este modelo está correcto y todos los estados son formalmente observables y controlables... Minreal no deberíae hacer nada... pero con la tolerancia por defecto ya hace algo porque este modelo es tan "poco" observable O controlable que  ya "dispara". Pero minreal no es la vía formalmente correcta de reducir orden en un modelo correctamente planteado; además, el resultado depende de tolerancias y cambios de variable en la realización (no es "equilibrada": balred reduce estados poco observables Y controlables).

size(minreal(barra40))
size(minreal(barra40,2e-11))

%% 2.1-Respuesta temporal/frecuencia del modelo de muelle (sin carga adicional)

step(barra40,1.5)
bodemag(barra40,logspace(0.5,3.5,250))
grid on

%% 2.2-Reducción basada en realización equilibrada del modelo de orden 80 (40 elementos)

[gg,baldata]=hsvd(barra40);
gg'
bar(gg)
semilogy(gg,'o'),grid on

% Nota: es tan poco controlable/observable que llegamos al límite de precisión numérica.
% Hagamos la reducción a orden 2 y 4 equilibrada.

reducido2=balred(barra40,2,baldata);
zpk(reducido2)
reducido4=balred(barra40,4,baldata);
zpk(reducido4)
pzmap(reducido4), grid on

% Comparemos también con el modelo muelle-amortiguador que descartara su masa y modos longitudinales, F=(K+fs)*x

s=tf('s');
Muelleamortiguador_sinmasa=1/(Ktotal+friccion_tot*s);
bodemag(barra40,reducido2,reducido4,Muelleamortiguador_sinmasa), grid on
legend('ord. 80','ord. 2','ord. 4','muelleamort sin masa')
err2=barra40-reducido2;
err4=barra40-reducido4;
bodemag(err2,err4), grid on
legend('ord. 2','ord. 4'), title('error de modelado (aditivo)')
norm(err2,Inf) %pico del error en frecuencia
norm(err4,Inf)
bodemag(err2/barra40,err4/barra40), grid on
legend('ord. 2','ord. 4'), title('error de modelado (relatlivo, porcentual)')

%% 2.3-Reducción basada en la separación de frecuencias (modal)

[Gslow,Gfast]=freqsep(barra40,100);
size(Gslow)
size(Gfast)
bodemag(barra40,Gslow,Gfast,logspace(0.5,3.5,200)), grid on
legend('completo','slow','fast')
bodemag(barra40,Gslow+dcgain(Gfast),Gfast-dcgain(Gfast),logspace(0.5,3.5,200)),grid on
legend('completo','slow+dcfit','fast-dcfit')

errslow=barra40-Gslow;
errslow_dc=barra40-Gslow-dcgain(Gfast);
bodemag(errslow,errslow_dc,err2,logspace(0.5,3.5,200)), grid on
legend('ord. 2modal','ord. 2modal+dcfit','balred'), title('error de modelado')
norm(errslow,Inf)
norm(errslow_dc,Inf)

%
%% 3.-Utilización del modelo reducido del muelle como subsistema
% Vamos a "cargar" el muelle con un objeto de masa M_2, y ejercer una fuerza F_2 sobra la masa grande M_2.
% Plantearemos las ecuaciones y veremos si todavía se puede reducir la complejidad del sistema más.
% Las ecuaciones combinadas serán: M_2\frac{d^2p}{dt^2}=F_2-F_{muelle}, junto con la dinámica entre la posición "p" y la fuerza del muelle anterior.
% En Trans. Laplace: M_2\cdot s^2\cdot p(s)=(F_2(s)-F_{muelle}(s)), p(s)=G_{muelle}(s)*F_{muelle}(s), con lo que resulta (M_2\cdot s^2+G^{-1}_{muelle}(s))\cdot p(s)=F_2(s)
% esto es, una FdT entre fuerza y posición dada por:
%  p(s)=\frac{1}{(M_2\cdot s^2+G^{-1}_{muelle}(s))}F_2(s)

M2=0.06;%probar con 0.06 [comparable a la del muelle] y con 1.2 Kg [mucho mayor que la del muelle].
ModeloTotalconmuellesinmasa=minreal(inv(M2*s^2+inv(Muelleamortiguador_sinmasa))) 
size(ss(ModeloTotalconmuellesinmasa))
ModeloTotalconmuelleorden2=minreal(inv(M2*s^2+inv(reducido2))); 
size(ModeloTotalconmuelleorden2)
ModeloTotalconmuelleorden4=minreal(inv(M2*s^2+inv(reducido4))); 
%aquí minreal sí está indicado, para eliminar errores de redondeo en manipulaciones
size(ModeloTotalconmuelleorden4)
ModeloTotalconmuelleorden80=minreal(inv(M2*s^2+inv(barra40))); %manejar sistemas de orden 82 complicado precisión numérica
size(ModeloTotalconmuelleorden80)
bodemag(ModeloTotalconmuelleorden2,ModeloTotalconmuelleorden4,ModeloTotalconmuelleorden80,ModeloTotalconmuellesinmasa,logspace(0,4,360)), grid on
legend('muelle2','muelle4','muelle80','muelle sin masa')

[gg2,baldata2]=hsvd(ModeloTotalconmuelleorden4);
gg2'
bar(gg2),grid on

% Parece que podemos eliminar dos de los cuatro estados... todo parece un masa-muelle de orden 2:

ModeloTotalRed=balred(ModeloTotalconmuelleorden4,2,baldata2);
zpk(ModeloTotalRed)
step(ModeloTotalconmuelleorden2,ModeloTotalconmuelleorden4,ModeloTotalRed,ModeloTotalconmuelleorden80,ModeloTotalconmuellesinmasa,6*M2)
legend('muelle2','muelle4','balred de muelle4','muelle80','muelle sin masa'), grid on
bodemag(ModeloTotalconmuelleorden4,ModeloTotalRed,ModeloTotalconmuelleorden80,ModeloTotalconmuellesinmasa,logspace(0,4,200)), grid on, 
legend('Orden 6 con din. muelle ord 4','Orden 2 (balred desde 6)','Con Din Muelle orden 80','ord 2 MuelleSinMasa','Location','best')

% Comparemos con el modelo más sencillo de no haber considerado NADA de dinámica del muelle:


bodemag(ModeloTotalconmuelleorden4,ModeloTotalRed,ModeloTotalconmuelleorden80,ModeloTotalconmuellesinmasa,logspace(1,2,250)), grid on
legend('Orden 6 con din. muelle ord 4','Orden 2 (balred desde 6)','Con Din Muelle orden 80','ord 2 MuelleSinMasa','Location','best')

%
%% Apéndice: funciones auxiliares

function sys=modelobarra(N,Ktotal, Mtotal, ftotal)
M=Mtotal/N;
K=Ktotal*N;f=ftotal*N;
%esto conserva las relaciones con fuerza de posición y velocidad suponiendo que cada elemento se mueve
%1/N veces lo que se mueve todo el muelle.

%estado 2N son [posiciones;velocidades];
A=zeros(N*2,N*2);
%entrada es fuerza extremo derecho (N)
B=zeros(N*2,1);
B(2*N)=1/M;
%salida es posición de la punta
C=zeros(1,N*2);
C(N)=1;
%deriv posicion son velocidades:
A(1:N,(N+1):(2*N))=eye(N);
%las aceleraciones tienen las ecuaciones siguientes:
for i=1:N
    if(i==1) %elemento pegado al extremo fijo
        A(N+i,i)=-2*K/M;
        A(N+i,N+i)=-2*f/M;
        A(N+i,i+1)=K/M;
        A(N+i,N+i+1)=f/M;
    elseif(i<N) %elementos intermedios
        A(N+i,i)=-2*K/M;
        A(N+i,N+i)=-2*f/M;
        A(N+i,i+1)=K/M;
        A(N+i,N+i+1)=f/M;
        A(N+i,i-1)=K/M;
        A(N+i,N+i-1)=f/M;
    else %i==N, extremo donde se aplica fuerza
        A(N+i,i)=-K/M;
        A(N+i,N+i)=-f/M;
        A(N+i,i-1)=K/M;
        A(N+i,N+i-1)=f/M;
    end
end
sys=ss(A,B,C,0);
sys.InputName={'Fuerza'};
sys.OutputName={'DesplazamientoEnExtremo'};
end