马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
这是一个优化实验 本实验有10个状态(1.5e-5 0 0 0 0 0.8 0 0.9 0 0)
11个参数(100000 1000 100 10000 50000 200 1000 20000 5000 100 500)，
我已得到初始系统的的远fisher矩阵（s'*s）， 现在要将其中一个状态的初始值(0.9)优化, 即想要让该状态在0.09到9之间变化并求出对应的Fisher矩阵。 问题是现在的运行结果不正确 请高手帮忙！
Fisher的程序代码：
clear all
clc;
format long;
for KN=1:100
IKK_normal=0.9;% initional value of N
IKK=IKK_normal*KN/10;
Xinit=[1.5e-5 0 0 0 0 0.8 0 IKK 0 0]';
k=[100000 1000 100 10000 50000 200 1000 20000 5000 100 500]';
sys = mdlInit(@ODESenfun, k , Xinit, zeros(11,10));
Tspan = 0:60:6000;
smltn = smltnInit(@ode15s, Tspan);
n = length(Xinit);
m = length(k);
[t, X, S3] = mdlPlay(sys, smltn);
error=10;
disp('please wait....')
S2 = S2DParam(S3);
FIM=S2'*S2;
end
[vec, lam]=eig(FIM);
for j=1:11
lambda(j)=lam(j,j);
end
%############################optimal##############
oed_a(KN)=trace(inv(FIM));
oed_d(KN)=det(FIM)
oed_e(KN)=min(lambda);
oed_me(KN)=max(lambda)/min(lambda);
问题就在上面的程序 希望高手帮忙下！！！
另外附上我得到初始fisher矩阵的程序！ 以下程序已经验证过 肯定正确！
mdllint.m:
function mdl = mdlInit(ode, theta, x_init, S_init)
% MDLINIT Initialise the parameters used in the model
% Validate input/output arguments
if nargin~=4
error('Function requires 4 input arguments');
end
% Initialise model data structure
mdl.ode = ode;
mdl.theta = theta;
mdl.x_init = x_init;
mdl.S_init = S_init;
mdl.thetaAeq = [];
mdl.thetaBeq = [];
mdlPlay.m:
function [t, X, S] = mdlPlay(mdl, smltn)
% MDLPLAY Simulate the model over the time span
%
% Returns the simulated states and the sensitivity matrix
%
% [t, X, S] = mdlPlay(mdl, smltn)
% mdl model data structure
% smltn simulation data structure
% X 2D matrix i,j (states, time) NB this is the transpose of normal
% return from the solver
% S 3D matrix i,j,k (states, parameters, time)
% Validate the input/output arguments
if nargin~=2,
error('Two input arguments are required');
end
if nargout<2 || nargout>3
error('Two or three output arguments are required');
end
% Initialise the parameters
num_states = length(mdl.x_init);
num_param = length(mdl.theta);
% Reshape the initial parameter/sensitivity values for the ODE simulation)
xs_init = [mdl.x_init; reshape(mdl.S_init, num_states*num_param, 1)];
% Run the ODE simulation
[t, X_S] = feval(smltn.solver, mdl.ode, smltn.tspan, xs_init, [], mdl.theta);
num_time = length(t);
% Reshape the state matrix (states, time)
X = X_S(:,1:num_states)';
% Reshape the sensitivity derivatives (states, parameters, time)
if nargout==3,
S = reshape(X_S(:,num_states+1:end)', num_states, num_param, num_time);
end
ODESenfun.m:
function dx = ODESenfun(t, x, k)
% ODESENFUN Solve both system dynamic ODEs and Sensitivity Equations
% System Dynamic ODEs
% x_dot = F*theta
% S_dot = J*S + F
% Check input and output arguments are correctly supplied
if nargin~=3
error('Three input arguments are required');
end
numStates = 10;
numPara = 11;
% Calculate the state derivatives
% Need to investigate a way to automatically generate these quantities (from a symbolic model)
%dx(1)=-k1*x(1)*x(6)+km1*x(2)+k4*x(4)-km4*x(1)*x(9)+k6*x(5);
%dx(2)=k1*x(1)*x(6) - km1*x(2) - k2*x(2) + km2*x(3)*x(7);
%dx(3)=k2*x(2) - km2*x(3)*x(7) - k3*x(3)*x(8) + km3*x(4) - k5W*x(3) + km5*x(5);
%dx(4)=k3*x(3)*x(8) - km3*x(4) - k4*x(4) + km4*x(1)*x(9);
%dx(5)=k5W*x(3) - km5*x(5) - k6*x(5);
%dx(6)=-k1*x(1)*x(6) + km1*x(2);
%dx(7)=k2*x(2) - km2*x(3)*x(7);
%dx(8)=-k3*x(3)*x(8) + km3*x(4);
%dx(9)=k4*x(4) - km4*x(1)*x(9);
%dx(10)=k6*x(5);
F = [-x(1)*x(6) x(2) 0 0 0 0 x(4) -x(1)*x(9) 0 0 x(5);
x(1)*x(6) -x(2) -x(2) x(3)*x(7) 0 0 0 0 0 0 0;
0 0 x(2) -x(3)*x(7) -x(3)*x(8) x(4) 0 0 -x(3) x(5) 0;
0 0 0 0 x(3)*x(8) -x(4) -x(4) x(1)*x(9) 0 0 0;
0 0 0 0 0 0 0 0 x(3) -x(5) -x(5);
-x(1)*x(6) x(2) 0 0 0 0 0 0 0 0 0;
0 0 x(2) -x(3)*x(7) 0 0 0 0 0 0 0;
0 0 0 0 -x(3)*x(8) x(4) 0 0 0 0 0;
0 0 0 0 0 0 x(4) -x(1)*x(9) 0 0 0;
0 0 0 0 0 0 0 0 0 0 x(5);];
dx = F*k;
% Calculate the sensitivity derivatives
% Need to investigate a way to automatically differentiate these quantities (from a symbolic model)
S = reshape(x(numStates+1:end), numStates, numPara);
J = [-k(1)*x(6)-k(8)*x(9) k(2) 0 k(7) k(11) -k(1)*x(1) 0 0 -k(8)*x(1) 0;
k(1)*x(6) -k(2)-k(3) k(4)*x(7) 0 0 k(1)*x(1) k(4)*x(3) 0 0 0;
0 k(3) -k(4)*x(7)-k(5)*x(8)-k(9) k(6) k(10) 0 -k(4)*x(3) -k(5)*x(3) 0 0;
k(8)*x(9) 0 k(5)*x(8) -k(6)-k(7) 0 0 0 k(5)*x(3) k(8)*x(1) 0;
0 0 k(9) 0 -k(10)-k(11) 0 0 0 0 0;
-k(1)*x(6) k(2) 0 0 0 -k(1)*x(1) 0 0 0 0;
0 k(3) -k(4)*x(7) 0 0 0 -k(4)*x(3) 0 0 0;
0 0 -k(5)*x(8) k(6) 0 0 0 -k(5)*x(3) 0 0;
-k(8)*x(9) 0 0 k(7) 0 0 0 0 -k(8)*x(1) 0;
0 0 0 0 k(11) 0 0 0 0 0;];
dS = J*S+F;
dx((numStates+1):(numStates*(numPara+1))) = reshape(dS,numStates*numPara,1);
S2DParam.m
function S1 = S2DParam(S)
% S2DPARAM Reshape the sensitivity matrix, S, to a 2D matrix
% S - 3D matrix (states, param, time)
% S1 - 2D matrix (states*time, param)
S1 = zeros(size(S,1)*size(S,3), size(S,2));
num_states = size(S,1);
num_time = size(S,3);
for i=1:num_states
S1((i-1)*num_time+1:i*num_time, :) = shiftdim(S(i, :, :),1)';
end
SensText.m
%% Calculate the 3-D sensitivity matrix (state-parameter-time),
% and transformed 2-D sensitivity matrix ((state*time)-parameter) for computing FIM
% This simulation is based on the ERK signal pathway model
clear all
clc
%%%% Set up initial model, true system and simulation parameters
Xinit=[1.5e-5 0 0 0 0 0.8 0 0.9 0 0]';
k = [100000 1000 100 10000 50000 200 1000 20000 5000 100 500]';
sys = mdlInit(@ODESenfun, k , Xinit, zeros(11,10));
% simulation time
Tspan = 0:1:6000;
smltn = smltnInit(@ode15s, Tspan);
n = length(Xinit);
m = length(k);
%%%% Solving the model states X and 3-D sensitivity matrix S
[t, X, S3] = mdlPlay(sys, smltn);
%%%% Transform the 3-D sensitivity matrix to 2-D sensitivity matrix
S2 = S2DParam(S3);
%S2(:,1) -- sensitivity coefficients for the 1st parameter. It include all
%the 10 states at all time points.
%S2(:,kth) -- sensitivity coefficients for the kth parameter. It include all
%the 10 states at all time points. (k=1,...,11)
%S2(1:101,kth) --- sensitivity coefficients of the 1st state with respect to the kth parameter
%S2(102:202,kth) --- sensitivity coefficients of the 2nd state with respect
%to the kth parameter
%tl=size(t);
%S2((n-1)*tl+1:n*tl,kth) --- sensitivity coefficients of the nth state with respect
%to the kth parameter
%n=1;
%m=11;
%s1=S2((n-1)*tl+1:n*tl,1:m)
%plot(t,s1)
tl=size(t);
con_str = ['\fontsize{20} I\kappaB\alpha (x1) ';...
'\fontsize{20} NF-\kappaB (x2) ';...
'\fontsize{20} I\kappaB\alpha-NF-\kappaB (x3) ';...
'\fontsize{20} I\kappaB\beta (x4) ';...
'\fontsize{20} I\kappaB\beta-NF-\kappaB (x5) ';...
'\fontsize{20} I\kappaB\epsilon (x6) ';...
'\fontsize{20} I\kappaB\epsilon-NF-\kappaB (x7) ';...
'\fontsize{20} IKKI\kappaB\alpha (x8) ';...
'\fontsize{20} IKKI\kappaB\alpha-NF-\kappaB (x9) ';...
'\fontsize{20} IKK (x10) '];
%%%%%%%%%%%figure out 10 state`s sensitivities with respect to 11 parameters
for n=1:10
figure(n)
s1=S2((n-1)*tl+1:n*tl,1);
s2=S2((n-1)*tl+1:n*tl,2);
s3=S2((n-1)*tl+1:n*tl,3);
s4=S2((n-1)*tl+1:n*tl,4);
s5=S2((n-1)*tl+1:n*tl,5);
s6=S2((n-1)*tl+1:n*tl,6);
s7=S2((n-1)*tl+1:n*tl,7);
s8=S2((n-1)*tl+1:n*tl,8);
s9=S2((n-1)*tl+1:n*tl,9);
s10=S2((n-1)*tl+1:n*tl,10);
s11=S2((n-1)*tl+1:n*tl,11);
plot(t,s1,t,s2,'y',t,s3,'b:',t,s4,'m-',t,s5,'r--',t,s6,'b-',t,s7,'g',...
t,s8,'m--',t,s9,'r-',t,s10,'k-.',t,s11,'k','Linewidth',2.5);
xlabel('\fontsize{16}time/s');
ylabel([con_str(n,:)]);
legend('k1','km1','k2','km2','k3','km3','k4','km4','k5W','km5','k6');
end
%%%%%%%%%%%%%%%%%%%%%%%%
FIM=S2'*S2;
smltnlnit.m
function smltn = smltnInit(solver, tspan)
% Initialize the parameters used in the ODE simulation
% Validate the input/output arguments
if nargin~=2
error('2 input arguments must be supplied');
end
% Set the data structure's fields
smltn.solver = solver;
smltn.tspan = tspan;
上面的所有程序中， 只有第一个程序错了 但是又不知道错在哪？ 请高手帮忙！ 鄙人不胜感激 ！！！
|
发表于 2011-3-22 00:53
提升卡
变色卡
显身卡
楼主
发表于 2011-3-22 23:58
