ADDED   mttroot/mtt/lib/control/PPP/ppp_RT.m
Index: mttroot/mtt/lib/control/PPP/ppp_RT.m
==================================================================
--- mttroot/mtt/lib/control/PPP/ppp_RT.m
+++ mttroot/mtt/lib/control/PPP/ppp_RT.m
@@ -0,0 +1,7 @@
+function [y,u] = ppp_RT (U)
+
+  ## usage:  [y,u] = ppp_RT (U)
+  ##
+  ## 
+
+endfunction

ADDED   mttroot/mtt/lib/control/PPP/ppp_RT_sim.m
Index: mttroot/mtt/lib/control/PPP/ppp_RT_sim.m
==================================================================
--- mttroot/mtt/lib/control/PPP/ppp_RT_sim.m
+++ mttroot/mtt/lib/control/PPP/ppp_RT_sim.m
@@ -0,0 +1,19 @@
+function [y,u] = ppp_RT_sim (U)
+
+  ## usage:  [y,u] = ppp_RT_sim (U)
+  ##
+  ## U   PPP weight (column vector)
+
+  global system_name_sim i_ppp_sim x_0_sim y_sim u_sim A_u_sim
+
+
+  ## Data from previous time - last point not used
+  if length(y_sim)>0		# Avoid initial junk
+    [n_t_old,junk] = size(y_sim);
+    y = y_sim(1:n_t_old-1,:); u = u_sim(1:n_t_old-1,:);
+  else
+    y=[]; u=[];
+  endif
+
+endfunction
+

ADDED   mttroot/mtt/lib/control/PPP/ppp_RT_sim_compute.m
Index: mttroot/mtt/lib/control/PPP/ppp_RT_sim_compute.m
==================================================================
--- mttroot/mtt/lib/control/PPP/ppp_RT_sim_compute.m
+++ mttroot/mtt/lib/control/PPP/ppp_RT_sim_compute.m
@@ -0,0 +1,26 @@
+function ppp_RT_sim_compute (U)
+
+  ## usage:  [y,u] = ppp_RT_sim_compute (U)
+  ##
+  ## U   PPP weight (column vector)
+
+  global system_name_sim i_ppp_sim x_0_sim y_sim u_sim A_u_sim simpar_sim
+
+  ## System details -- defines simulation within ol interval
+  par = eval(sprintf("%s_numpar;", system_name_sim));
+  t = [0:simpar_sim.dt:simpar_sim.last];
+  n_t = length(t);
+  [n_x,n_y,n_u] = eval(sprintf("%s_def;", system_name_sim));
+  [n_U,junk] = size(A_u_sim);
+
+  ## Set up u_star
+  u_star = ppp_ustar(A_u_sim,1,t,0,0,n_u-n_U);
+
+  ## Simulate
+  par(i_ppp_sim(:,3)) = U;		# Update the simulation ppp weights
+  [y_sim,x] = eval(sprintf("%s_sim(x_0_sim, par, simpar_sim, u_star);", \
+			     system_name_sim));
+  x_0_sim  = x(n_t,:)';     # Extract state for next time
+  u_sim = (u_star*U);
+endfunction
+

ADDED   mttroot/mtt/lib/control/PPP/ppp_nlin_run.m
Index: mttroot/mtt/lib/control/PPP/ppp_nlin_run.m
==================================================================
--- mttroot/mtt/lib/control/PPP/ppp_nlin_run.m
+++ mttroot/mtt/lib/control/PPP/ppp_nlin_run.m
@@ -0,0 +1,242 @@
+function [y,u,t,p,UU,t_open,t_ppp,t_est,its_ppp,its_est] = ppp_nlin_run (system_name,i_ppp,i_par,A_u,w_s,N_ol,extras)
+
+
+  ## usage: [y,u,t,p,U,t_open,t_ppp,t_est,its_ppp,its_est] =
+  ## ppp_nlin_run (system_name,i_ppp,i_par,A_u,w_s,N_ol[,extras])
+  ##
+  ##  y,u,t   System output, input and corresponding time p
+  ##  Estimated parameters U       PPP weight vector t_open  The
+  ##  open-loop interval t_ppp   Time for each ppp optimisation t_est
+  ##  Time for each estimation i_ppp   Matrix of ppp gain indices: col.
+  ##  1 Gain indices in sensitivity system col. 2  Gain sensitivity
+  ##  indices in sensitivity system col. 3  Gain indices in  system
+  ##  i_par Matrix of indices of estimated parameters col. 1  Parameter
+  ##  indices in sensitivity system col. 2  Parameter sensitivity
+  ##  indices in sensitivity system col. 3  Parameter indices in  system
+  ##  A_u     Basis function generating matrix w_s     w_star: That part
+  ##  of the moving horizon setpoint within the optimisation horizon.
+  ##  N_ol The number of open-loop intervals to be computed extras
+  ##  Extra parameters in  a structure: extras.alpha           ??
+  ##  extras.criterion Optimisation convergence criterion
+  ##  extras.emulate_timing Simulate some real-time features
+  ##  extras.estimate Estimate parameters and states
+  ##  extras.max_iterations Maximum optimisation iterations
+  ##  extras.simulate 1 for simulation (not real-time) extras.vInitial Levenberg-Marquardt parameter
+  ##          extras.verbose         1 for extra info display
+  ##
+  ## Real-time implementatipn of  nonlinear PPP Copyright (C) 2001,2002
+  ## by Peter J. Gawthrop
+
+  ## Globals to pass details to simulation
+  global system_name_sim i_ppp_sim x_0_sim y_sim u_sim A_u_sim simpar_sim
+
+  ## Defaults
+  if nargin<7
+    extras.alpha = 0.1;
+    extras.criterion = 1e-5;
+    extras.emulate_timing = 0;
+    extras.estimate = 1;
+    extras.max_iterations = 10;
+    extras.simulate = 1;
+    extras.v = 1e-5;
+    extras.verbose = 0;
+  endif
+  
+  ## Names
+  s_system_name = sprintf("s%s", system_name);
+
+  ## System details -- defines simulation within ol interval
+  par = eval(sprintf("%s_numpar;", system_name));
+  simpar = eval(sprintf("%s_simpar;", system_name));
+  dt = simpar.dt;		# Sample interval
+  simpar_est = simpar;		# Initial estimation simulation params
+  simpar_model = simpar;	# Initial internal model simulation params
+  simpar_pred = simpar;		# Initial prediction simulation params
+  T_ol_0 = simpar.last;		# The initial specified interval
+  n_t = round(simpar.last/simpar.dt); # Corresponding length
+  x_0 = eval(sprintf("%s_state(par);", system_name));
+  x_0_model = x_0;
+  [n_x,n_y,n_u] = eval(sprintf("%s_def;", system_name));
+
+  ## Sensitivity system details -- defines moving horizon simulation
+  simpars = eval(sprintf("%s_simpar;", s_system_name));
+  pars = eval(sprintf("%s_numpar;", s_system_name));
+  x_0s = eval(sprintf("%s_state(pars);", s_system_name));
+
+  ## Times
+  ## -- within opt horizon
+  n_Tau = round(simpars.last/simpars.dt);
+  dtau = simpars.dt;
+  Tau = [0:n_Tau-1]'*dtau;
+  [n_tau,n_w] = size(w_s);
+  tau = Tau(n_Tau-n_tau+1:n_Tau);
+  w = w_s(n_tau,:);		# Final value of setpoint
+
+
+  ## Main simulation loop
+  y = [];
+  x = [];
+  u = [];
+  t = [];
+
+  p = [];
+
+  t_last = 0;
+  UU = [];
+  UU_l =[];
+  UU_c =[];
+  
+  t_ppp = [];
+  t_est = [];
+  its_ppp = [];
+  its_est = [];
+  t_open = [];
+  x_nexts = zeros(2*n_x,1);
+
+  ## Initial U is zero
+  [n_U,junk] = size(A_u);
+  U = zeros(n_U,1);
+
+  ## Create input basis functions
+  u_star_tau = ppp_ustar(A_u,1,Tau',0,0,n_u-n_U);
+	
+  ## Reverse time to get "previous" U
+  U_old = U;
+
+  if (extras.simulate==1)
+    ## Set up globals for simulation
+    system_name_sim = system_name;
+    i_ppp_sim = i_ppp;
+    x_0_sim = x_0;
+    y_sim = [];			# Junk
+    u_sim = [];			# Junk
+    A_u_sim = A_u;
+    simpar_sim = simpar;
+    T_total = simpar.last;
+  endif
+  
+  for i = 0:N_ol		# Main loop 
+    printf("%i",i);
+    if (extras.simulate==1)
+      [y_ol,u_ol] = ppp_RT_sim(U); # Simulate
+    else
+      [y_ol,u_ol] = ppp_RT(U);	# Real thing
+    endif
+
+    t_start = time;		# Initialise time
+
+    if (i==0)			# Data is rubbish at i=0 - ignore
+      usleep(T_ol_0*1e6);		# Hang about
+    else
+      ## Set up time information for the gathered data
+      n_ol = length(y_ol);	# How many data points?
+      t_ol = [0:n_ol-1]'*dt;
+      T_ol = n_ol*dt;		# Length of ol interval
+      t_open = [t_open;T_ol];
+
+      ## Generate input to actual system
+      u_star_t = ppp_ustar(A_u,1,t_ol',0,0,n_u-n_U);
+
+      ## Tune parameters/states
+      if (extras.estimate==1)
+	## Save up the estimated parameters
+	par_est = pars(i_par(:,1));
+	p = [p; par_est'];
+
+	## Set up according to interval length
+	if (T_ol>T_ol_0) ## Truncate data
+	  simpar_est.last = T_ol_0;
+	  y_est = y_ol(1:n_t+1,:);
+	else
+	  simpar_est.last = T_ol;
+	  y_est = y_ol;
+	endif
+
+	simpar_pred.last = T_ol_0; # Predicted length of next interval
+	pars(i_ppp(:,1)) = U_old; # Update the simulation ppp weights
+	
+	## Optimise
+	tick = time;
+	[pars,Par,Error,Y,its] = \
+	    ppp_optimise(s_system_name,x_0s,pars,simpar_est,u_star_t,y_est,i_par,extras);
+	est_time = time-tick;  
+	t_est = [t_est;est_time];
+	its_est = [its_est; its-1];
+      endif
+
+      ## Update internal model
+      par(i_ppp(:,3)) = U_old; # Update the simulation ppp weights
+
+      if (extras.estimate==1)
+	par(i_par(:,3)) = pars(i_par(:,1)); # Update the simulation params
+      endif
+      
+      simpar_model.last = T_ol;
+      [y_model,x_model] = eval(sprintf("%s_sim(x_0_model, par, simpar_model, \
+ 					       u_star_t);",system_name));
+
+      x_0 = x_model(n_ol+1,:)';	# Initial state of next interval
+      x_0_model = x_0;
+
+      ## Compute U by optimisation
+      tick = time;
+
+      ## Predict state at start of next interval
+      par(i_ppp(:,3)) = U;
+      [y_next,x_next] = eval(sprintf("%s_sim(x_0, par, simpar, \
+					     u_star_t);",system_name));
+      x_next = x_next(n_t+1,:)'; # Initial state for next time
+      x_nexts(1:2:(2*n_x)-1) = x_next;
+      
+      ## Optimize for next interval      
+      U_old = U;		# Save previous value
+      U = expm(A_u*T_ol)*U;	# Initialise from continuation trajectory
+      pars(i_ppp(:,1)) = U;	# Put initial value of U into the parameter vector
+      [U, U_all, Error, Y, its] = ppp_nlin(system_name,x_nexts,pars,simpars,u_star_tau,w_s,i_ppp,extras);
+
+
+      ppp_time = time-tick;  
+      t_ppp = [t_ppp;ppp_time];
+      its_ppp = [its_ppp; its-1];
+
+      ## Total execution time
+      T_total = time - t_start;
+      if (extras.simulate==1)&&(extras.emulate_timing!=1)
+	printf(".");
+	T_diff = 0;		# Always correct interval
+      else
+	T_diff = T_total - T_ol_0; # Compute difference
+	if T_diff<0
+	  printf("-");
+	  usleep(-T_diff*1e6);
+	  T_total = time - t_start;
+	else
+	  printf("+");
+	endif   
+				printf("%2.2f",T_total);
+      endif
+      T_total = simpar.dt*round(T_total/simpar.dt); # Whole no. of intervals
+
+      pars(i_ppp(:,1)) = U;	# Put final value of U into the parameter vector
+
+      ## Save up data
+      y_ol = y_ol(1:n_ol,:);
+      y = [y; y_ol];
+      u = [u; u_ol];
+      UU = [UU; U'];
+      t = [t; t_ol+t_last*ones(n_ol,1) ];
+      t_last = t_last + T_ol;
+
+    endif
+
+
+    if (extras.simulate==1) # Do the actual simulation
+      if (extras.emulate_timing==1) # Emulate timing
+	simpar_sim.last = T_total; # simulate for actual execution time
+      endif
+      ppp_RT_sim_compute (U_old);
+    endif
+    
+  endfor
+  printf("\n");
+endfunction