%% Identificación experimental de modelo 1er-orden+integrador a partir de datos de BUCLE CERRADO (identificación "indirecta")
% © 2025, Antonio Sala Piqueras, Universitat Politècnica de València, Spain. Todos los derechos reservados.
% Presentaciones en vídeo:
% https://personales.upv.es/asala/YT/V/idbcind1.html  [intro, motivación, planteamiento]
% https://personales.upv.es/asala/YT/V/idbcind2.html  [ejemplo escalón referencia]
% https://personales.upv.es/asala/YT/V/idbcind3.html [ejemplo escalón variable manipulada]
%
% Objetivos: a veces se pueden hacer experimentos sólo con un controlador funcionando, bien porque no se desee "desconectar" un proceso "en producción", bien porque el proceso sea inestable y la única forma de operarlo sea con un controlador activo. Aquí se presentan un par de ejemplos sencillos de cómo podría abordarse el problema de hacer un ensayo de identificación del proceso, suponiendo que el controlador K(s) es conocido. Esto se denomina "identificación en bucle cerrado". Lo haremos con identificación "indirecta": identificando una función de transferencia de bucle cerrado y luego "despejando" la del proceso G en bucle abierto.
%
% NOTA: lo que aquí contamos son ejemplos "de juguete" (el "toy example" en inglés) para un primer curso de grado sobre control automático... la identificación en bucle cerrado es más compleja que la de "bucle abierto", en general, ver nota/conclusion al final.
% Tabla de Contenidos
% Modelo del proceso
% Ejemplo 1: identificación a partir de bucle con control P ante escalón en referencia
%     Generación de experimento simulado
%     Identificación experimental 2o orden indirecta
%     Recuperación de FdT del proceso en bucle abierto G
% Ejemplo 2: identificación ante ensayo de "incremento de variable manipulada" con control PD.
%     Generación de experimento simulado
%     Identificación experimental 2o orden indirecta
%     Recuperación de la FdT del proceso
%% Modelo del proceso
%
% *Supondremos T_carga=0.
% \frac{d\omega}{dt}=\frac{1}{J}(T_{mot}-b\cdot \omega), o sea Js\cdot \omega(s)=T_{mot}(s)-b\cdot \omega(s), esto es: \omega(s)=\frac{1}{Js+b}T_{mot}(s).
% \frac{d\theta}{dt}=\omega, o sea s\cdot \theta(s)=\omega(s), esto es:
% \theta(s)=\frac{1}{s}\omega(s)=\frac{1}{s\cdot (Js+b)}T_{mot}(s).

s=tf('s');
J=0.1;b=0.11;
Gvel=1/(J*s+b);
Gpos=Gvel/s;
zpk(Gpos)

% ID bucle cerrado: se desea identificar G(s) con datos de la salida y(s) de este bucle de control, con K(s) conocido:
%
%
%
%
%
% \bullet La identificación a partir del par (u(t),y(t)) se denominará identificación "directa"... Tiene ciertos problemas si hay ruidos ( \delta_u desconocida, ruido de medida)... porque u entonces no es "accesible", sólo tenemos u_k; entonces,  con r=0 pero perturbaciones actuando, lo que se tiene es u_K=-Ky, o sea que identificamos y=\frac{-1}{K}u_K.
% \bullet La identificación a partir del par (r(t),y(t)) o del par (\delta_u(t),y(t)) se denominará identificación "indirecta". Identificamos una función de bucle cerrado que es, a su vez, función de G y K... luego, como K es conocido, despejamos G. Ese es el objetivo de este material.
%% Ejemplo 1: identificación a partir de bucle con control P ante escalón en referencia
%
%% Generación de experimento simulado
% La salida ante escalón en referencia de un bucle de control estándar es:
% y(s)=\frac{G(s)\cdot K(s)}{1+G(s)\cdot K(s)}\cdot r(s).
% En efecto, y=G\cdot u=GK\cdot e=GK\cdot (r-y), con lo que (1+GK)\cdot y=GK\cdot r, y tenemos la fórmula buscada.

K=0.25;%P    
%K=0.25+0.05*s; %PD sale con "ceros" en T, con lo que no valen las fórmulas
%de abajo... Habría que mirar fórmulas más complicadas o estimar "por
%mínimos cuadrados" como "procest" o similar.

T=feedback(Gpos*K,1);
zpk(T)

% Estos serían los únicos datos disponibles

h1=stepplot(T,Gpos*K,9); grid on
h1.Responses(1).LineWidth=2;ylim([0 2])
legend("Bucle cerrado P","bucle abierto")
title("Respuesta a escalón unitario de REFERENCIA")

%% Identificación experimental 2o orden indirecta
% Se trata de un sistema de 2o orden oscilatorio, sin ceros, por lo que aplicaremos las fórmulas que dicen que el sistema
% G(s)=\frac{[Gain]\cdot \omega_n^2}{s^2+2\xi\omega_n\cdot s+\omega_n^2 , en caso 0<\xi<1 tiene una respuesta ante escalón unitario (c.i. nulas) dada por:
%
% con A=1, siendo \sigma=\xi\omega_n (parte real polos), \omega_p=\omega_n\sqrt{1-\xi^2} (parte imaginaria polos).
% Su valor máximo es y_{max}=(1+\delta)K, con \delta=e^{\frac{-\xi\pi}{\sqrt{1-\xi^2}}, y el primer pico se alcanza en \frac{\pi}{\omega_p} segundos. El tiempo de establecimiento 98% es \approx\frac{4}{\sigma}=\frac{4}{\xi\omega_n}.
%
% Redibujemos y marquemos datos "experimentales".

hb=stepplot(T); grid on
hb.Responses.LineWidth=2;

% Experimentalmente, medimos

y_final=1; y_inicial=0;
ymax=1.31;
yline(ymax,'m',Label="Y_{max}");
t_pico=2.1; 
xline(t_pico,'r',Label="1er pico",LabelVerticalAlignment="bottom")

% Deberíamos hacer una primera comprobación que t_pico es el semiperíodo de las oscilaciones, o sea que el período asociado a \omega_p es 4.2 segundos.

xline(3*t_pico,'r',Label="1er pico+1 período",LabelVerticalAlignment="bottom")

% Con lo que, una vez verificado, aplicaríamos las fórmulas:

delta=(ymax-y_final)/(y_final-y_inicial)
tmp=log(delta)^2;
xi=sqrt(tmp/(tmp+pi^2)) %coeficiente de amortiguamiento
wp=pi/t_pico %frecuencia propia
wn=wp/sqrt(1-xi^2)
Gain=y_final/1 %incremento salida / incremento entrada, en caso general
%NO puedo llamarle 'K' porque ya está la letra usada para el controlador.

% La identificación estaría terminada, pero comprobaremos que el tiempo de establecimiento observado cuadra con los cálculos antes de dar todo por concluido.

sigma=xi*wn %parte real polos
t_est=4/sigma %comprobar para confirmar que cálculos correctos
xline(t_est,'g',Label="t_{est}",LabelVerticalAlignment="bottom")

% Todo cuadra, con lo que tenemos nuestro estimado:

T_estim=Gain*wn^2/(s^2+2*xi*wn*s+wn^2)

%% Recuperación de FdT del proceso en bucle abierto G

T=T_estim

% De T=GK/(1+GK) tenemos T+GKT=GK, T=GK(1-T), G=\frac{T}{K\cdot (1-T)}=\frac{1}{K\cdot (T^{-1}-1)}

Gestim=1/K/(1/T-1);
zpk(Gestim)

% ¡Lo tenemos!
%% Ejemplo 2: identificación ante ensayo de "incremento de variable manipulada" con control PD.
% A veces el ensayo es con "escalón en entrada al proceso" aparte del control... Podría hasta incluso ser más intuitivo: doy la entrada justo a lo que quiero identificar.
%% Generación de experimento simulado
% En ese caso, y(s)=\frac{G(s)}{1+G(s)K(s)}\cdot \delta_u(s); en efecto, y=G\cdot (\delta_u+K(r-y)), o sea, si r no tiene incrementos, y=G\cdot (\delta_u-Ky), con lo que (1+GK)y=G\delta_u, y tenemos la expresión buscada.
% En bucle abierto hubiese sido y(s)=G(s)\cdot \delta_u(s), evidentemente.

K=0.25+0.06*s; %PD
W=feedback(Gpos,K)
h=stepplot(W); grid on
h.Responses.LineWidth=2;
title("Respuesta a escalón unitario en \delta_u")

%% Identificación experimental 2o orden indirecta
% Las mismas fórmulas (2o orden subamortiguado) que en caso anterior siguen aplicando.
% Repetimos plot para marcar datos experimentales

h2=stepplot(W,9);grid on
h2.Responses.LineWidth=2;

% Experimentalmente, medimos

y_final=4; y_inicial=0;
ymax=4.55;
yline(ymax,'m',Label="Y_{max}");
t_pico=2.35; 
xline(t_pico,'r',Label="1er pico",LabelVerticalAlignment="bottom")

% Deberíamos hacer una primera comprobación que t_pico es el semiperíodo de las oscilaciones, o sea que el período asociado a \omega_p es 4.2 segundos.

xline(3*t_pico,'r',Label="1er pico+1 período",LabelVerticalAlignment="bottom")

% Con lo que, una vez verificado, aplicaríamos las fórmulas:

delta=(ymax-y_final)/(y_final-y_inicial)
tmp=log(delta)^2;
xi=sqrt(tmp/(tmp+pi^2)) %coeficiente de amortiguamiento
wp=pi/t_pico %frecuencia propia
wn=wp/sqrt(1-xi^2)
Gain=(y_final-y_inicial)/1 %incremento salida / incremento entrada, en caso general
%NO puedo llamarle 'K' porque ya está la letra usada para el controlador.

% La identificación estaría terminada, pero comprobaremos que el tiempo de establecimiento observado cuadra con los cálculos antes de dar todo por concluido.

sigma=xi*wn %parte real polos
t_est=4/sigma %comprobar para confirmar que cálculos correctos
xline(t_est,'g',Label="t_est",LabelVerticalAlignment="bottom")

% Todo cuadra, con lo que tenemos nuestro estimado:

W_estim=Gain*wn^2/(s^2+2*xi*wn*s+wn^2)

%% Recuperación de la FdT del proceso
% y(s)=\frac{G(s)}{1+G(s)K(s)}\cdot \delta_u(s). Si denomino y(s)=W(s)\cdot \delta_u(s), entonces \frac{G}{1+GK}=W resulta en G=W+GKW, esto es (1-KW)G=W, o sea G=\frac{W}{1-KW}=\frac{1}{W^{-1}-K}.

Gestim2=minreal(1/(1/W_estim-K))

% Nota: todo está muy `preparado', y esto debe considerarse únicamente un "ejercicio académico". Dependiendo de la estructura del regulador y del modelo del proceso, pueden salir otros tipos de funciones de transferencia de orden superior o con ceros, y las saturaciones, otras nolinealidades, ruidos de proceso y de medida, etc. pueden hacer que este procedimiento de resultados bastante malos en algunas situaciones reales: la identificación en bucle cerrado no es un problema `trivial' en la práctica.