帮助将代码从MATLAB转换为python

Question

我正在将一些代码从MATLAB翻译成python。这段代码模拟了模型的行为，我想从中估计参数。问题是，使用python和MATLAB获得的结果非常不同。我决定使用GEKKO，而不是使用scipy，但我是个新手，所以我发现了一些问题。我报告了MATLAB代码和scipy代码。 我报告了我的代码的最小示例

ef cost_function(parameters, measured, t, N, total_active):
    Tp = simulate(parameters, t, N, total_active,measured[0])
    measured = np.squeeze(measured)
    Tp = Tp[:,2]
    return (np.sum(np.power((np.divide(np.diff(measured), measured[1:])-np.divide(np.diff(Tp),Tp[1:])),2)))


def SIR_model(x, t, params, total_active):
    S0, I0, R0 = x
    v, tau, _ = params
    dSdt = - v * S0 * I0 / total_active(t)
    dIdt = v * S0 * I0 / total_active(t) - g * I0 - tau * I0
    dCdt = tau * I0
    return [dSdt, dIdt, dCdt]


def simulate(p, t, N, total_active, Q0):
    T = np.zeros((len(t), 3))
    x0 = [N - 100 * p[2] - Q0, 100*p[2], Q0]
    T[0, :] = [N - 100 * p[2] - Q0, 100*p[2], Q0]

    for i in range(len(t) - 1):
        ts = [t[i], t[i + 1]]
        y = odeint(SIR_model, x0, ts, args=(p, total_active))
        x0 = y[-1]
        T[i + 1, :] = x0
    return T

def identify_model(data, initial_guess, t, N, total_active):
    # Set bounds
    It0 = initial_guess['It0']/100
    v = initial_guess['v']
    tau = initial_guess['tau']
    lim_sup = [v, tau * 1.5, (It0 * 1.3)]
    lim_inf = [v * 0, tau * 0.9, (It0 * 0.7)]
    #Set initial condition of the model
    bounds = Bounds(lim_inf, lim_sup)
    # estimate parameters in order to minimize cost function
    options = {"maxiter": 100, "ftol": 1e-06}  # with SLSQP
    cons = {'type': 'ineq', 'fun': constraints, "args": (N, total_active, t,data[0])}
    return minimize(cost_function, np.asarray([v, tau, 100*It0]), args=(data, t, N, total_active), bounds=bounds,
                    constraints=cons, options=options, tol=1e-06,
                    method="SLSQP")['x']

data=[275.5,317.,457.33333333,646.,888.66666667,1236.66666667,1619.33333333,2077.33333333,2542.33333333]
times = [i for i in range(len(data))]
total_active_data=[59999725.,59999684.33333334,59999558.66666666,59999385.33333334,59999158.33333333,59998823.,59998474.66666666,59998053.33333333,59997652.66666666]
total_active = interp1d([i for i in range(len(total_active_data))], total_active_data, fill_value="extrapolate")
initial_guess = {"v": 0.97, "tau": 0.066, "It0": 100}
print(identify_model(data,initial_guess,times,60e6,total_active))

结果，这个片段给出了[0.97099097,0.099,130.]。 （我希望如此）等效的MATLAB代码是：

function [pars] = Identify_Model(data,initial_guess,lim_inf,lim_sup,times,N,total_active)

scalefac=100;

%Initialize the initial guess for each parameter
v=initial_guess.v;
It0=initial_guess.It0;
tau=initial_guess.tau;
g=0.07;

lim_inf(3)=lim_inf(3)/scalefac;
lim_sup(3)=lim_sup(3)/scalefac;

%Identify the model parameters
options=optimset('MaxIter',100,'TolFun',1e-6,'TolX',1e-6);
parmin=fmincon(@(pars) error_SIR(data,[pars(1),pars(2),g,scalefac*pars(3)],times,N,total_active),[v,tau,It0],[],[],[],[],lim_inf,lim_sup,[],options);
pars=num2cell([parmin(1:2),scalefac*parmin(3)]);

pars=[pars(1),pars(2),g,pars(3)];

end



function [costo]=error_SIR(data,pars,tempi,N,totale_attivi)

%Assign the parameters
pars=num2cell(pars);
[v,tau,g,I0]=pars{:};
%Simulate the model
Q0=data(1,1);
S0=N-I0-Q0;

[~,y]=ode45(@(t,x) SIR(t,x,v,tau,g,totale_attivi),tempi,[S0;I0;Q0]);
if length(tempi)==2
    y=[y(1,:);y(end,:)];
end

%Identify on the normalized data (Data/mean data)
costo=sum(((diff(sum(data,1))./sum(data(:,2:end),1))-(diff(sum(y(:,3),2))./sum(y(2:end,3),2))').^2);
end

function y=SIR(t,x,v,tau,g,total_active)
    y=zeros(3,1);
    y(1)=-v*x(1)*x(2)/total_active(t);           % S
    y(2)=v*x(1)*x(2)/total_active(t)-(tau+g)*x(2);   % I
    y(3)=tau*x(2);      % C
end

total_active_data=[59999725.,59999684.333,59999558.666,59999385.333,59999158.33,59998823.,59998474.66,59998053.333,59997652.66666666]
total_active  = @(t) interp1(1:9,total_active_data,t); 
initial_guess.It0=100;
initial_guess.v=0.97;
initial_guess.g=0.07;  
initial_guess.tau=0.066;
g=initial_guess.g;
It0=initial_guess.It0;
v=initial_guess.v;
tau=initial_guess.tau;
N=60e6
%Define the constrints for the identification
lim_sup=[v*10, tau*1.5, It0*1.3];
lim_inf=[v*0, tau*0.9, It0*0.7];
data=[275.5,317.,457.33333333,646.,888.66666667,1236.66666667,1619.33333333,2077.33333333,2542.33333333]

times=1:9;
%identify the model parameters
pars=Identify_Model(data,initial_guess,lim_inf,lim_sup,times,N,total_active)

如果我想用GEKKO翻译，我面临几个问题：

• 如何在GEKKO中设置成本函数
• 如何获得预测结果
• 由于两个参数可以随时间变化，我开发了一种统计方法来了解它是否发生了变化，我可以直接计算MHE吗？怎么做

我已经用GEKKO编写了这段代码，希望它是正确的：

m = GEKKO()
m.time = [i for i in range(4)]
t=[i for i in range(4)]
v = m.FV(0.97, 0, 9.7)
tau = m.FV(0.066, 0.0594, 0.099)
I0 = m.FV(100, 70, 130)
I0.STATUS = 1
tau.STATUS = 1
v.STATUS = 1
S0 = m.CV(60e6 - 275.5 - 100)
I0 = m.CV(100)
R0= m.CV(275.5)
total_active_data=[59999725.,59999684.33333334,59999558.66666666,59999385.33333334]
total_active_data = m.Param([59999725.,59999684.33333334,59999558.66666666,59999385.33333334])
m.Equation(S0.dt()==-v*I0*S0/total_active_data)
m.Equation(I0.dt() == v * S0 * I0/total_active_data - 0.07 * I0 - tau * I0)
m.Equation(R0.dt()==0.07*I0)
m.options.IMODE = 5  # MHE
# Solve
m.solve(disp=True)

Answer 1

你有一些很接近的东西。您可以用m.Minimize()定义目标，例如用平方误差和目标匹配数据

S_data = m.Param([60e6,59.999e6,59.888e6,59.777e6])
m.Minimize((S0-S_data)**2)

要计算I0的初始条件，请使用函数m.free_initial(I0)。其他FV定义正确。我为拟合创建了一些样本数据（蓝色圆点）

from gekko import GEKKO
m = GEKKO(remote=False)
m.time = [i for i in range(4)]
t=[i for i in range(4)]
v = m.FV(0.97, lb=0, ub=9.7); v.STATUS = 1
tau = m.FV(0.066, lb=0.0594, ub=0.099); tau.STATUS = 1

S0 = m.CV(60e6 - 275.5 - 100)
I0 = m.CV(100,lb=0)
m.free_initial(I0)
R0= m.CV(275.5)

total_active_data=[59999725.,59999684.33333334,\
                   59999558.66666666,59999385.33333334]
total_active_data = m.Param([59999725.,59999684.33333334,\
                             59999558.66666666,59999385.33333334])
S_data = m.Param([60e6,59.999e6,59.888e6,59.777e6])
m.Equation(S0.dt()==-v*I0*S0/total_active_data)
m.Equation(I0.dt()== v*S0 * I0/total_active_data - 0.07*I0 - tau*I0)
m.Equation(R0.dt()== 0.07*I0)
m.options.IMODE = 5  # MHE

# Objective
m.Minimize((S0-S_data)**2)

# Solve
m.solve(disp=True)

import matplotlib.pyplot as plt
plt.subplot(3,1,1)
plt.ylabel('S')
plt.plot(m.time,S0.value,'b ')
plt.plot(m.time,S_data.value,'bo')
plt.subplot(3,1,2)
plt.ylabel('I')
plt.plot(m.time,I0.value,'r ')
plt.plot([0],[100],'ro',markersize=10)
plt.subplot(3,1,3)
plt.ylabel('R')
plt.plot(m.time,R0.value,'k:')
plt.xlabel('time')
plt.show()