lphtoby 发表于 2011-3-22 00:53

急!!matlab程序上的问题, 请高手帮忙!

这是一个优化实验 本实验有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=';
    k=';
    sys = mdlInit(@ODESenfun, k , Xinit, zeros(11,10));
    Tspan = 0:60:6000;
    smltn = smltnInit(@ode15s, Tspan);
    n = length(Xinit);
    m = length(k);
    = mdlPlay(sys, smltn);
    error=10;
    disp('please wait....')
    S2 = S2DParam(S3);
    FIM=S2'*S2;
end

=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 = mdlPlay(mdl, smltn)
% MDLPLAY        Simulate the model over the time span
%
% Returns the simulated states and the sensitivity matrix
%
% = 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 = ;
% Run the ODE simulation
= 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)
% S2DPARAMReshape 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=';
k = ';

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
= 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 11parameters
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();
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;




上面的所有程序中, 只有第一个程序错了 但是又不知道错在哪? 请高手帮忙!鄙人不胜感激 !!!

lphtoby 发表于 2011-3-22 00:59

回复 1 # lphtoby 的帖子

这是10个状态11个变量的系统! 怎么得到其中一个状态变在一定范围内变化化对应的fisher矩阵
高人也可以就本方面给予指导谢谢!!

ChaChing 发表于 2011-3-22 23:58

提问的智慧!!!!(发帖前请认真阅读)
http://forum.vibunion.com/forum/viewthread.php?tid=21991
:@)
页: [1]
查看完整版本: 急!!matlab程序上的问题, 请高手帮忙!