function [par,Par,Error,Y,iterations] = ppp_optimise(system_name,x_0,u,t,par_0,free,y_0,extras);
## Usage: [par,Par,Error,Y,iterations] = ppp_optimise(system_name,x_0,u,t,par_0,free,y_0,weight,extras);
## system_name String containg system name
## y_s actual system output
## theta_0 initial parameter estimate
## free Indices of the free parameters within theta_0
## weight Weighting function - same dimension as y_s
## extras.criterion convergence criterion
## extras.max_iterations limit to number of iterations
## extras.v Initial Levenberg-Marquardt parameter
## Copyright (C) 1999,2000 by Peter J. Gawthrop
## 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_y,n_data] = size(y_0);
n_th = length(i_s);
error_old = inf;
error_old_old = inf;
error = 1e50;
reduction = 1e50;
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
while (abs(reduction)>extras.criterion)&&\
(abs(error)>extras.criterion)&&\
(iterations<extras.max_iterations)
iterations = iterations + 1; # Increment iteration counter
[y,y_par] = eval(sprintf("%s_sim(x_0,u,t,par,i_s)", system_name)); # Simulate
##Evaluate error, cost derivative J and cost second derivative JJ
error = 0;
J = zeros(n_th,1);
JJ = zeros(n_th,n_th);
for i = 1:n_y
E = y(i,:) - y_0(i,:); # Error
error = error + (E*E'); # Sum the error over outputs
y_par_i = y_par(i:n_y:n_y*n_th,:);
J = J + y_par_i*E'; # Jacobian
JJ = JJ + y_par_i*y_par_i'; # Newton Euler approx Hessian (with LM
# term
endfor
if iterations>1 # 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 extras.verbose
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
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
endwhile
endfunction