function [par,Par,Error,Y,iterations] = \ ppp_optimise(system_name,x_0,par_0,simpar,u,y_0,free,extras); ## Levenberg-Marquardt optimisation for PPP/MTT ## Usage: [par,Par,Error,Y,iterations] = ppp_optimise(system_name,x_0,par_0,simpar,u,y_0,free[,extras]); ## system_name String containing system name ## x_0 Initial state ## par_0 Initial parameter vector estimate ## simpar Simulation parameters: ## .first first time ## .dt time increment ## .stepfactor Euler integration step factor ## u System input (column vector, each row is u') ## y_0 Desired model output ## free one row for each adjustable parameter ## first column parameter indices ## second column corresponding sensitivity indices ## extras (opt) optimisation parameters ## .criterion convergence criterion ## .max_iterations limit to number of iterations ## .v Initial Levenberg-Marquardt parameter ###################################### ##### Model Transformation Tools ##### ###################################### ############################################################### ## Version control history ############################################################### ## $Id$ ## $Log$ ## Revision 1.2 2001/04/04 08:36:25 gawthrop ## Restuctured to be more logical. ## Data is now in columns to be compatible with MTT. ## ## Revision 1.1 2000/12/28 11:58:07 peterg ## Put under CVS ## ############################################################### ## Copyright (C) 1999,2000 by Peter J. Gawthrop sim_command = sprintf("%s_ssim(x_0,par,simpar,u,i_s)", system_name); ## Extract indices i_t = free(:,1); # Parameters i_s = free(:,2)'; # Sensitivities if nargin<8 extras.criterion = 1e-5; extras.max_iterations = 10; extras.v = 1e-5; extras.verbose = 0; endif [n_data,n_y] = size(y_0); n_th = length(i_s); error_old = inf; error_old_old = inf; error = 1e50; reduction = inf; predicted_reduction = 0; par = par_0; Par = par_0; step = ones(n_th,1); Error = []; Y = []; iterations = 0; v = extras.v; # Levenverg-Marquardt parameter. r = 1; # Step ratio if extras.verbose # Diagnostics printf("Iteration: %i\n", iterations); printf(" error: %g\n", error); printf(" reduction: %g\n", reduction); printf(" prediction: %g\n", predicted_reduction); printf(" ratio: %g\n", r); printf(" L-M param: %g\n", v); printf(" parameters: "); for i_th=1:n_th printf("%g ", par(i_t(i_th))); endfor printf("\n"); endif while (abs(reduction)>extras.criterion)&&\ (abs(error)>extras.criterion)&&\ (iterations1 # Adjust the Levenberg-Marquardt parameter reduction = error_old-error; predicted_reduction = 2*J'*step + step'*JJ*step; r = predicted_reduction/reduction; if (r<0.25)||(reduction<0) v = 4*v; elseif r>0.75 v = v/2; endif if reduction<0 # Its getting worse par(i_t) = par(i_t) + step; # rewind parameter error = error_old; # rewind error error_old = error_old_old; # rewind old error if extras.verbose printf(" Rewinding ....\n"); endif endif endif ## Compute step using pseudo inverse JJL = JJ + v*eye(n_th); # Levenberg-Marquardt term step = pinv(JJL)*J; # Step size par(i_t) = par(i_t) - step; # Increment parameters error_old_old = error_old; # Save old error error_old = error; # Save error ##Some diagnostics Error = [Error error]; # Save error Par = [Par par]; # Save parameters Y = [Y y]; # Save output if extras.verbose # Diagnostics printf("Iteration: %i\n", iterations); printf(" error: %g\n", error); printf(" reduction: %g\n", reduction); printf(" prediction: %g\n", predicted_reduction); printf(" ratio: %g\n", r); printf(" L-M param: %g\n", v); printf(" parameters: "); for i_th=1:n_th printf("%g ", par(i_t(i_th))); endfor printf("\n"); endif endwhile endfunction