Cómo ajustar un conjunto de datos a una función definida en Simulink

Tenemos un conjunto de puntos. Por ejemplo, sean los puntos:

tm = [1 2 3 4 6 8 10 15]';
ym = [4.8 3.1 2.2 1.3 0.5 0.2 0.1 0.0]';

Y queremos aproximarlos a la siguiente función:

y' = -k·y

y(0) = y0

Para lo que tenemos que estimar los parámetros "k" e "y0".

Como la función contiene una derivada, lo más cómodo es definiarla usando Simulink. Para ello, creamos el siguiente archivo "simu1.mdl" en Simulink:

- El bloque "Fcn", donde hemos definido la función, está en "User-Defined Functions".

- El bloque "Integrator" está  en "Continuous". Donde digo que la "Initial condition source" es "external".

- Los bloques "Constant" y "Clock" están en "Sources", mientras que el bloque "To Workspace" está en "Sinks".

- En el bloque "To Workspace" defino el formato de salida como "Array" (no "Structure" que es lo que viene por defecto). Además, defino el nombre de la variables de salida como "yout".

- En tiempo de simulación hemos puesto la variable "tend". Antes de poder ejecutar el programa, habrá que definir el valor de "tend", "k" e "y0" en Matlab.

- Cuando ejecutamos este programa, Matlab, además de calcular el vector "yout" con el valor de la variable "y", también calcula el vector "tout", con el valor del tiempo.

---

A continuación creamos la función objetivo a minimizar, que calcula la suma de los errores al cuadrado, para unos valores de "k" y "y0" dados. A este archivo le damos el nombre "fobj1.m".

function z = fobj01(tm,ym,u)
global k y0

k = u(1);
y0 = u(2);

sim('simu1.mdl');
ycal = resample(timeseries(yout,tout),tm);
ycal = ycal.data;
clc

z = sum((ym - ycal).^2);

---

Finalmente, desde Matlab, ya podemos hacer la estimación de parámetros. Para ello escribimos:

tm = [1 2 3 4 6 8 10 15]';
ym = [4.8 3.1 2.2 1.3 0.5 0.2 0.1 0.0]';
tend = 20;

global k y0
[u,f] = fminsearch(@(u) fobj01(tm,ym,u),[1 1])

donde "u" es el valor de los parámetros estimados, y "f" es el valor de función objetivo (suma de los error al cuadro.)

Matlab nos devuelve los siguientes valores:

u =

    0.4247    7.3344


f =

    0.0309

---

Por último, podemos dibujar tanto los valores medidos como los estimados:

k = u(1);
y0 = u(2);
sim('simu1.mdl');

figure(1)
hold on
box on
plot(tout,yout,'-','Color',[1 0.6 0.78],'LineWidth',2)
plot(tm,ym,'o','MarkerSize',5,'Color',[0.48 0.06 0.89],'MarkerFaceColor',[0.48 0.06 0.89])
lg1 = legend({'real','estimado'},1,'EdgeColor',[1 1 1]);
axis([0 20 0 8])
set(gca,'XTick',[0:5:20])

set(gca,'YTick',[0:2:8])

set(gca,'Fontsize',10)