function [y,x,u,t,U,U_c,U_l] = ppp_nlin_sim (system_name,A_u,tau,t_ol,N_ol,w,extras)

  ## usage:  [y,x,u,t,U,U_c,U_l] = ppp_nlin_sim (system_name,A_u,tau,t_ol,N,w)
  ##
  ## 
  
  ## Simulate nonlinear PPP
  ## Copyright (C) 2000 by Peter J. Gawthrop

  ## Defaults
  if nargin<7
    extras.U_initial = "zero";
    extras.U_next = "continuation";
    extras.criterion = 1e-5;
    extras.max_iterations = 10;
    extras.v = 0.1;
    extras.verbose = 0;
  endif
  
  

  ## Names
  s_system_name = sprintf("s%s", system_name);

  ## System details 
  par = eval(sprintf("%s_numpar;", system_name));
  x_0 = eval(sprintf("%s_state;", system_name))
  [n_x,n_y,n_u] = eval(sprintf("%s_def;", system_name));

  ## Sensitivity system details
  sympars = eval(sprintf("%s_sympar;", s_system_name));
  pars = eval(sprintf("%s_numpar;", s_system_name));

  ## Times
  n_tau = length(tau);
  dtau = tau(2)-tau(1);		# Optimisation sample time
  Tau = [0:dtau:tau(n_tau)+eps]; # Time  in the moving axes
  n_Tau = length(Tau);
  dt = t_ol(2)-t_ol(1);
  n_t = length(t_ol);
  T_ol = t_ol(n_t)+dt
  

  ## Weight and moving-horizon setpoint
  weight = [zeros(n_y,n_Tau-n_tau), ones(n_y,n_tau)]; 
  ws = w*ones(1,n_tau);

  ## Create input basis functions
  [n_U,junk] = size(A_u);

  ## For moving horizon
  eA = expm(A_u*dtau);
  u_star_i = ones(n_U,1);
  u_star_tau = [];
  for i = 1:n_Tau
    u_star_tau = [u_star_tau, u_star_i];
    u_star_i = eA*u_star_i;
  endfor

  ## and for actual implementation
  eA = expm(A_u*dt);
  u_star_i = ones(n_U,1);
  u_star_t = [];
  for i = 1:n_t
    u_star_t = [u_star_t, u_star_i];
    u_star_i = eA*u_star_i;
  endfor
  
  if extras.verbose
    title("U*(tau)")
    xlabel("tau");
    plot(Tau,u_star_tau)
  endif
  

  ## Check number of inputs adjust if necessary
  if n_u>n_U
    disp(sprintf("Augmenting inputs with %i zeros", n_u-n_U));
    u_star_tau = [u_star_tau; zeros(n_u-n_U, n_Tau)];
    u_star_t   = [u_star_t; zeros(n_u-n_U, n_t)];
  endif
  
  if n_u<n_U
    error(sprintf("n_U (%i) must be less than or equal to n_u (%i)", \
		  n_U, n_u));
  endif
		  
  ## Indices of U and sensitivities
  FREE = "[";
  free_name = "ppp";
  for i=1:n_U
    FREE = sprintf("%ssympars.%s_%i sympars.%s_%is", FREE, free_name, \
		   i, free_name, i);
    if i<n_U
      symbol = "; ";
    else
      symbol = "];";
    endif
    FREE = sprintf("%s%s", FREE, symbol);
  endfor
FREE
  free = eval(FREE);

  ## Compute linear gains 
  [A,B,C,D] = eval(sprintf("%s_sm(par);", system_name));
  B = B(:,1); D = D(:,1);
  [k_x,k_w,K_x,K_w] = ppp_lin(A,B,C,D,A_u,0,tau);

  ## Main simulation loop
  y = [];
  x = [];
  u = [];
  t = [];
  t_last = 0;
  UU = [];
  UU_l =[];
  UU_c =[];
  
  x_0s = zeros(2*n_x,1);

  if  strcmp(extras.U_initial,"linear")
    U = K_w*w - K_x*x_0;
  elseif strcmp(extras.U_initial,"zero")
    U = zeros(n_U,1);
  else
    error(sprintf("extras.U_initial = \"%s\" is invalid", extras.U_initial));
  endif

  ## Reverse time to get "previous" U
   U = expm(-A_u*T_ol)*U;

  for i = 1:N_ol
    ## Compute initial U from linear case
    U_l = K_w*w - K_x*x_0;

    ## Compute initial U  for next time from continuation trajectory
    U_c = expm(A_u*T_ol)*U;

    ## Create sensitivity system state
    x_0s(1:2:2*n_x) = x_0;
    
    ## Set next U (initial value for optimisation)
    if  strcmp(extras.U_next,"linear")
      U = U_l;
    elseif strcmp(extras.U_next,"continuation")
      U = U_c;
    elseif strcmp(extras.U_next,"zero")
      U = zeros(n_U,1);
    else
      error(sprintf("extras.U_next = \"%s\" is invalid", extras.U_next));
    endif
    ## Put initial value of U into the parameter vector
    pars(free(:,1)) = U;

    ## Compute U
    tick = time;
    if extras.max_iterations>0
      [U, U_all, Error, Y] = ppp_nlin(s_system_name,x_0s,u_star_tau,tau,pars,free,ws,extras);
    else
      Error = [];
    endif
    opt_time = time-tick;  
    printf("Optimisation %i took %i iterations and %2.2f sec\n", i, \
	   length(Error), opt_time);
#     title(sprintf("Moving horizon trajectories: Interval %i",i));
#     grid;
#     plot(tau,Y)
  
    ## Generate control

    u_ol = U'*u_star_t(1:n_U,:);

    ## Simulate system 
    ## (Assumes that first gain parameter is one)
    U_ol = [u_ol; zeros(n_u-1,n_t)]; # Generate the vector input
    [y_ol,x_ol] = eval(sprintf("%s_sim(x_0, U_ol, t_ol, par);", system_name));

    ## Extract state for next time
    x_0  = x_ol(:,n_t);

    y = [y y_ol(:,1:n_t)];
    x = [x x_ol(:,1:n_t)];
    u = [u u_ol(:,1:n_t)];

    UU = [UU, U];
    UU_l = [UU_l, U_l];
    UU_c = [UU_c, U_c];

    t = [t, t_ol(:,1:n_t)+t_last*ones(1,n_t) ];
    t_last = t_last + t_ol(n_t); 
  endfor

  ## Rename returned matrices
  U = UU;
  U_l = UU_l;
  U_c = UU_c;
endfunction