"initial_state" had to be defined.
globals had to be included where the variables were needed.
"y" has been removed from the globals because it seems to be a solution variable.
"u" has to be used as "u(t,x)", not simply as "u" in function "modelo_stolwijk_hardy".
"main" has to be called from the script part.
ode45 has been replaced by ode15s because your ODE system seems to be stiff.
"grid on" has been added to the first subplot.
I don't remember if maybe I forgot some more changes.
% Parámetros del modelo de Stolwijk y Hardy
global k1 k2 k3 c theta alpha betaa
k1 = 0.003;
k2 = 0.006;
k3 = 0.001;
c = 100;
theta = 2.5;
alpha = 7600;
betaa = 139000;
main()
function main()
% Parámetros del controlador PID
Kp = 0.1;
Ki = 0.05;
Kd = 0.01;
int_error = 0;
prev_error = 0;
% Tiempo de simulación
tspan = [0 180]; % Desde 0 a 180 minutos
% Función de referencia (nivel deseado de glucosa)
reference = @(t) 100;
% Función de bomba de insulina
u = @(t, x) bomba_insulina(t, x);
initial_state = [1;1];
% Simulación del sistema
[t, x] = ode15s(@(t, x) modelo_stolwijk_hardy(t, x, u), tspan, initial_state);
% Resultados
glucosa = x(:, 1);
insulina = x(:, 2);
% Gráficos
clf();
figure;
subplot(211);
plot(t, glucosa);
xlabel('Tiempo (minutos)');
ylabel('Glucosa en sangre (mg/dL)');
title('Concentración de Glucosa en Sangre');
grid on
subplot(212);
plot(t, insulina);
xlabel('Tiempo (minutos)');
ylabel('Insulina en sangre (mU/L)');
title('Concentración de Insulina en Sangre');
% Mostrar las gráficas
grid on;
end
function dx = modelo_stolwijk_hardy(t, x, u)
global k1 k2 k3 c theta alpha betaa
G = x(1);
y = x(2);
% Utilizar la variable t
dx = [k1 * (u(t,x) - G) - k2 * G;
betaa * (G - theta) - alpha * y];
end
function u = bomba_insulina(t, x)
global k1 k2 k3 c theta alpha betaa
G = x(1);
y = x(2);
% Condición de no producción de insulina
if G >= theta
u = 0;
else
u = betaa * (G - theta) - alpha * y;
end
end
