function [theta,Theta,Error,Y,iterations] = mtt_optimise (system_name,y_s,theta_0,method,free,weight,criterion,max_iterations,alpha)
  ## Usage: [theta,Theta,Error,Y,iterations] = mtt_optimise (system_name,y_s,theta_0,method,free,weight,criterion,max_iterations,alpha)
  ##  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
  ##  method "time" or "freq"
  ##  criterion convergence criterion
  ##  max_iterations limit to number of iterations
  ##  alpha  Optimisation gain parameter

  ## Copyright (C) 1999 by Peter J. Gawthrop

  if nargin<4
    method="time";
  endif

  N = length(theta_0);
  if nargin<5
    free = [1:N];
  endif
  
  if nargin<6
    weight = ones(size(y_s));
  endif

  if nargin<7
    criterion = 1e-5;
  endif
  
  if nargin <8
    max_iterations = 10;
  endif

  if nargin<9
    alpha = 1.0;
  endif
  
  if (!strcmp(method,"time"))&&(!strcmp(method,"freq"))
    error("method must be either time or freq")
  endif
  
  N_theta = length(free);
  Weight = weight*ones(1,N_theta); # Sensitivity weight
  e_last = 1e20;
  error=1e10;
  theta = theta_0;
  Theta = [];
  Error = [];
  Y = [];
  iterations = -1;
  while (abs(e_last-error)>criterion)&&(iterations<max_iterations)
    iterations = iterations + 1;
    e_last = error;
    eval(sprintf("[t,y,y_theta] = mtt_s%s(system_name,theta,free);",method)); # Simulate system
    Theta = [Theta theta];	# Save parameters
    Y = [Y y];			# Save output
    E = weight.*(y - y_s);	# Weighted error
    y_theta = Weight.*y_theta;	# Weighted sensitivity
    error = (E'*E);		# Sum the error
    Error = [Error error];
    ##    theta(free) = theta(free) - alpha*(real(y_theta'*y_theta)\real(y_theta'*E));
    tol = 1e-4;
    JJ = real(y_theta'*y_theta);
    ## sigma = svd(JJ)
    theta(free) = theta(free) - alpha*( pinv(JJ,tol)*real(y_theta'*E) );
  endwhile

endfunction