%% Comprensión "intuitiva" del comportamiento de un sistema dinámico lineal: cálculo de tren de escalones para comportamiento deseado
%
% © 2022, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Presentación en vídeo en http://personales.upv.es/asala/YT/V/linregla3.html
% Este código ejecutó sin errores en Matlab R2022b (Linux)
% Objetivo: comprender intuitivamente el concepto de sistema (dinámico) lineal invariante en el tiempo, y cómo usar la respuesta temporal para calcular perfiles de entrada aproximados que consigan ciertos objetivos.
% Tabla de Contenidos
% Ensayo escalón del que obtener datos
% Resolución de problemas de cálculo de escalones de entrada
%     Valor final prefijado
%     Valor final y tiempo de establecimiento prefijado
% Apéndice funciones auxiliares
%
%% Ensayo escalón del que obtener datos
% Consideremos un cierto sistema desconocido, del que sólo sabemos que es lineal (o aproximadamente) alrededor de un punto de operación u=2, y=4;
% Simulemos su respuesta ante escalón de 0.5 (incremental)

u_pf=2; y_pf=4;
inc_u=0.5;
u=@(t) u_pf+inc_u*(t>=0);
Y=simulasistema(u);
yline(y_pf,'-.g',Label="Punto de funcionamiento de la salida")
yline(u_pf,'-.c',Label="Punto de funcionamiento de la entrada")

% Marcamos en la gráfica ciertas lineas que necesitaremos en los ejemplos posteriores

xline(1.4,'-.') %lo necesitamos luego
yline(5.25,'-.r') %también lo necesitamos luego
yline(4.66,'-.') %también lo necesitamos luego
legend("Salida","Entrada"), title("Experimento 1: escalón")
valfin=Y(end) %valor final

%% Resolución de problemas de cálculo de escalones de entrada
% Con sólo esa información del ensayo ante escalón, podemos contestar a varias preguntas.
%% Valor final prefijado
% 1.) ¿Qué entrada será necesaria para subir a 5.25 unidades la salida?
% Por "proporcionalidad", si con incremento entrada 0.5 sube 0.75 unidades, para subir a 5.25:

incsalida=valfin-y_pf
gananciaestatica=incsalida/inc_u %incremento por unidad de entrada
incsalida_deseado=5.25-y_pf
inc_entrada_calculado=incsalida_deseado/gananciaestatica

% Por tanto, el valor "absoluto" que en equilibrio conseguirá la salida deseada será:

u_calculado=u_pf+inc_entrada_calculado
u=@(t) u_pf+inc_entrada_calculado*(t>=0);
simulasistema(u);
xline(1.4,'-.')
yline(5.25,'-.r',Label="valor final deseado",LabelHorizontalAlignment="left")
legend("Salida","Entrada",Location="best")

%% Valor final y tiempo de establecimiento prefijado
% 1.) ¿Qué entrada será necesaria para subir a 5.25 unidades la salida en 1.4 segundos?
% Por "proporcionalidad", si con incremento 0.5 sube 0.66 unidades en 1.4 segundos, para subir a 5.25:

incsalida=4.66-y_pf
ganancia_en_1punto4segundos=incsalida/inc_u %incremento por unidad de entrada
incsalida_deseado=5.25-y_pf
inc_entrada_calculado2=incsalida_deseado/ganancia_en_1punto4segundos
u=@(t) u_pf+inc_entrada_calculado2*(t>=0);
simulasistema(u);
yline(5.25,'-.r',Label="valor deseado",LabelHorizontalAlignment="left")
xline(1.4,'-.',Label="tiempo de establecimiento deseado",LabelVerticalAlignment="bottom")
legend("Salida","Entrada",Location="best")

% Para que no "se pase" habrá que conmutar al valor de entrada que mantiene el valor deseado en equilibrio una vez se alcanza:

u=@(t) u_pf+inc_entrada_calculado2*(t>=0).*(t<=1.4)+inc_entrada_calculado*(t>1.4);
simulasistema(u);
yline(5.25,'-.r',Label="valor deseado",LabelHorizontalAlignment="left")
xline(1.4,'-.',Label="tiempo de establecimiento deseado",LabelVerticalAlignment="middle")
legend("Salida","Entrada",Location="best")

% Esto es, hemos diseñado un perfil de 2 escalones: un "boost" inicial para que suba más rápido y luego un valor "permanente" para mantenerse en el punto deseado.
% NOTA: los cálculos realizados sólo son "exactos" en sistemas lineales de primer orden; en sistemas lineales de orden superior hay cierta "inercia" que hará que aunque se baje la entrada hasta el punto de equilibrio calculado, haya un cierto "sobrepasamiento transitorio". Esto se ilustrará en un segundo vídeo http://personales.upv.es/asala/YT/V/linregla3ord2.html
%
%
%
%% Apéndice funciones auxiliares
% Este código es "secreto": no hace falta examinarlo para realizar los cálculos objetivo de este material.

function dxdt=modelo1(x,u)
 dxdt=-1.5*x+2.25*u+1.5;
end

function Y=simulasistema(u)
opts=odeset(Reltol=1e-5,AbsTol=1e-5);
[T,X]=ode45(@(t,x) modelo1(x,u(t)),[-1 8],4,opts);
plot(T,X,LineWidth=2), grid on
Y=X;
hold on
plot(T,u(T'),'--',LineWidth=2,Color=[.3 .75 .65])
hold off
end

%