#!/bin/csh
## Automatically generated from bashrc on Wed May  9 09:02:08 BST 2001 - DO NOT EDIT
#! /bin/sh

     ###################################### 
     ##### Model Transformation Tools #####
     ######################################

# Bourne shell script: mttrc - sets up paths etc for mtt
# Usage: mttrc 

# Copyright (c) P.J.Gawthrop 1996,1977.


###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.22  2001/04/12 03:08:00  geraint
## Improved sh->csh conversion, reduces environment namespace pollution.
##
## Revision 1.21  2001/04/10 13:56:13  gawthrop
## Uses standard mkoctfile
##
## Revision 1.20  2001/04/10 13:08:19  gawthrop
## Smoother translation to .cs using sh2csh
##
## Revision 1.19  2001/03/30 15:13:49  gawthrop
## Rationalised simulation modes to each return mtt_data
##
## Revision 1.18  2001/03/19 02:28:52  geraint
## Branch merge: merging-ode2odes-exe back to MAIN.
##
## Revision 1.17.2.4  2001/03/06 03:48:43  geraint

## Initial merge of ode2odes.exe into main mtt.
## standalone_rep.make deleted: rules moved to mtt, variables to mttrc.
##
## Revision 1.17  2000/12/27 16:46:13  peterg
## Stripped the mtt- from paths
##
## Revision 1.16  2000/12/27 15:16:44  peterg

##
## Revision 1.15  2000/12/27 14:57:43  peterg
## Now takes the base path as an argument
##
## Revision 1.14  2000/12/27 13:11:43  peterg
## *** empty log message ***
##

## Revision 1.7  1998/03/24 09:11:49  peterg
## Compatible with .csh version
##
## Revision 1.6  1998/03/13 11:53:29  peterg
## reduce --> reduce 64
##
## Revision 1.5  1998/01/16 08:55:01  peterg
## MAKE=make
##
## Revision 1.4  1998/01/06 09:14:51  peterg
## Added latex2html to setup
##
# Revision 1.3  1998/01/06  09:11:26  peterg
# Removed matlab from the setup
#
# Revision 1.2  1997/12/04  10:49:16  peterg
# Put under RCS at last
# Added CC variable
#
###############################################################

## When using csh, replace /home/peterg/Development/mttroot/mtt by the mtt base path, eg /usr/share/mtt/latest
setenv MTT_BASE /home/peterg/Development/mttroot/mtt




  echo Setting paths with base $MTT_BASE
  # The following line sets up the make to use -- gmake is the standard 
  # but you may wish to use lsmake for parallelism
  setenv MAKE 'make'
  
  # The following sets up the c compiler
  setenv CC 'gcc'

  # Setup latex2html
  setenv LATEX2HTML "latex2html -contents_in_navigation -index_in_navigation -address http://mtt.sourceforge.net"
  
  # Ascend stuff
  setenv ASCENDLIBRARY $MTTPATH/ascend/lib
  
  # Oct file generation - use version with no optimisation.
  #setenv MKOCTFILE $MTT_LIB/octave/mkoctfile # This for no optimisation
    setenv MKOCTFILE mkoctfile

  # ode2odes.exe stuff

    # local system

set PLAT="i686-pc-linux-gnu"
#    PREFIX="/usr/local"
set PREFIX="/usr"
set GCCVERS="2.95.2"
set SRCOCTAVE="/cvs/octave"

#    PLAT="mips-sgi-irix6.5"
#    PREFIX="/usr/people/bevangp/GNU"
#    GCCVERS="2.95.2"
#    SRCOCTAVE="${PREFIX}/../build/octave-2.1.33"

    # include paths

set IOCTAVE="-I${PREFIX}/include/octave"

    # library paths

#    LOCTAVE="-L${PREFIX}/lib/octave -loctave -lcruft -loctinterp"
set LOCTAVE="-L${PREFIX}/lib/octave -loctave -lcruft -loctinterp"
set LKPATHSEA="-L${SRCOCTAVE}/kpathsea -lkpathsea"
set LREADLINE=" -L${SRCOCTAVE}/readline -lreadline"
set LSYSTEM="-ldl -lm -lncurses"
set LF2C="-L${PREFIX}/lib/gcc-lib/${PLAT}/${GCCVERS} -lg2c"

    # compiler options

set DEBUG="-g"
set OPTIM="-O3"
set FLAGS="-fno-rtti -fno-exceptions -fno-implicit-templates"

    # exported variables

    setenv MTT_CXX "g++"
    setenv MTT_CXXFLAGS "${DEBUG} ${OPTIM} ${FLAGS}"
    setenv MTT_CXXLIBS "${LOCTAVE} ${LKPATHSEA} ${LREADLINE} ${LF2C} ${LSYSTEM}"
    setenv MTT_CXXINCS "-I. ${IOCTAVE}"
    setenv MTT_LDFLAGS " "

echo Creating $outfile
make_m() {
mtt_header ${sys} ssim m >  ${outfile}
cat >> ${outfile} <<EOF

  ## Pass input as a global
  global MTT_input MTT_input_index MTT_input_last

  if nargin>4  
    [n,m] = size(index);
    if (n>m)       # Make sure its a row vector
      index = index';
    endif
  endif

  [nx2,ny2] = ${sys}_def;

  y_par = [];
  MTT_input = u;
  [MTT_input_last,m] = size(u);

  if nargin<5			# No index given
    MTT_input_index = 0;
    [mtt_data] = ${sys}_ode2odes(x0,par,simpar);
    y = mtt_data(:,2:2:1+ny2);
    x = mtt_data(:,3+ny2:2:2+ny2+nx2);
    ypar = [];
  else				# Compute sensitivities as well
    for i=index
      MTT_input_index = 0;
      p = par;           # Reset parameters
      p(i) = 1;          # Set sensitivity index to 1
      [mtt_data] = ${sys}_ode2odes(x0,p,simpar);
      if (i==index(1))
	y = mtt_data(:,2:2:1+ny2);
      endif
      y_par = [y_par, mtt_data(:,3:2:2+ny2)];
      x = [];
    endfor
  endif

endfunction
EOF

}

# Copyright (C) 2000 by Peter J. Gawthrop

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.32  2001/05/24 07:42:12  gawthrop
## Included and updated the missing tf_r2m
##
## Revision 1.31  2001/04/03 14:49:42  gawthrop
## Revised to incorporate new ssim (sensitivity simulation)
## representation (m only just now).
##
## Revision 1.30  2001/03/30 15:13:58  gawthrop
## Rationalised simulation modes to each return mtt_data
##

        args=mttpar
        declarestates=yes
	;;
    ssim)
	states=no;
        inputs=no;
	parameters=no;
        output='y,y_par,x'
        args='x0,par,simpar,u,index'
	;;
    tf)
	states=no;
	inputs=no;
	parameters=yes;
        output='mttnum,mttden'

%% Label file for system sCS (sCS_lbl.txt)
%SUMMARY sCS C component with initial state


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.1  2000/12/28 10:31:35  peterg
% %% Put under RCS
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Port aliases
%ALIAS	in	in_1,in_2

% Argument aliases
%ALIAS	$1	effort,c
%ALIAS	$2	e_0
%ALIAS	$3	c_s
%ALIAS	$4	e_0s

%% Each line should be of one of the following forms:
%	     a comment (ie starting with %)
%	     component-name	cr_name	arg1,arg2,..argn
%	     blank

% ---- Component labels ----

% Component type C
	c		slin	effort,c;c_s

% Component type SS
	[in]	 SS	external,external

% Component type Se
	e_0		SS	e_0;e_0s

function [U, U_all, Error, Y] = ppp_nlin (system_name,x_0,par_0,sim,us,w,free,extras)

  ## usage:  [U, U_all, Error, Y] = ppp_nlin (system_name,x_0,par_0,sim,us,w,free,extras)
  ##
  ## 

   if nargin<8
     extras.criterion = 1e-8;
     extras.max_iterations = 10;
     extras.v = 0.1;
     extras.verbose = 1;
   endif

   s_system_name = sprintf("s%s", system_name); # Name of sensitivity system

   ## Details
   [n_x,n_y,n_u] = eval(sprintf("%s_def;", system_name));
   [n_tau,n_us] = size(us);

   ## Checks
   if (n_us<>n_u)
     error(sprintf("Inputs (%i) not equal to system inputs (%i)", n_us, n_u));
   endif
   
  [par,Par,Error,Y] = ppp_optimise(s_system_name,x_0,par_0,sim,us,w,free,extras);

  U = par(free(:,1));
  U_all = Par(free(:,1),:);
endfunction

function [y,x,u,t,UU,UU_c,UU_l] = ppp_nlin_sim (system_name,i_ppp,i_par,A_u,w_s,N_ol,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

  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))
  x_0 = eval(sprintf("%s_state(par);", system_name))
  [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))
  sympars = eval(sprintf("%s_sympar;", s_system_name));
  pars = eval(sprintf("%s_numpar;", 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(length(w_s));		# Final value of setpoint


  ## -- within ol interval
  n_t = round(simpar.last/simpar.dt)
  dt = simpar.dt;
  t_ol = [0:n_t-1]'*dt;
  T_ol = n_t*dt;

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

  ## 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)
    title("U*(t)")
    xlabel("t_ol");
    plot(t_ol,u_star_t)
  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
  
















  ## 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;

  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		# Main loop 
    ## 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
    pars(i_ppp(:,1)) = U;	# Put initial value of U into the parameter vector


    ## Compute U by optimisation
    tick = time;
    if extras.max_iterations>0
      [U, U_all, Error, Y] = ppp_nlin(system_name,x_0s,pars,simpars,u_star_tau,w_s,i_ppp,extras);
      pars(i_ppp(:,1)) = U;	# Put final value of U into the parameter vector
    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_star_t*U;		# Not used - just for show


    ## Simulate system over one ol interval
    ## (Assumes that first gain parameter is one)
    ##U = [u_ol zeros(n_t,n_u-1)]; # Generate the vector input
    [y_ol,y_s,x_ol] = eval(sprintf("%s_ssim(x_0s, pars, simpar, u_star_t);", s_system_name));
    
    x_0  = x_ol(n_t+1,:)';	# Extract state for next time
    y_ol = y_ol(1:n_t,:);	# Avoid extra points due to rounding error 
    x_ol = x_ol(1:n_t,:);	# Avoid extra points due to rounding error 


    y = [y; y_ol];
    x = [x; x_ol];
    u = [u; u_ol];

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

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





endfunction

  ######################################
  
  ###############################################################
  ## Version control history
  ###############################################################
  ## $Id$
  ## $Log$
  ## Revision 1.3  2001/04/05 11:50:12  gawthrop
  ## Tidied up documentation + verbose mode
  ##
  ## 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
  ##

    extras.max_iterations = 10;
    extras.v = 1e-5;
    extras.verbose = 0;
  endif
  

  [n_data,n_y] = size(y_0);
  if n_data<n_y
    error("ppp_optimise: y_0 should be in columns, not rows")
  endif

  n_th = length(i_s);
  error_old = inf;
  error_old_old = inf;
  error = 1e50;
  reduction = inf;
  predicted_reduction = 0;

  while (abs(reduction)>extras.criterion)&&\
	(abs(error)>extras.criterion)&&\
	(iterations<extras.max_iterations)

    iterations = iterations + 1; # Increment iteration counter

    [y,y_par] = eval(sim_command); # Simulate
    [N_data,N_y] = size(y);

    if (N_y!=n_y)
      mess = sprintf("n_y (%i) in data not same as n_y (%i) in model", n_y,N_y);
      error(mess);
    endif
    
    ## Use the last part of the simulation to compare with data
    y = y(1+N_data-n_data:N_data,:);
    y_par = y_par(1+N_data-n_data:N_data,:);

    ##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

function Ustar = ppp_ustar (A_u,n_u,tau,order,packed)

  ## usage:  Us = ppp_ustar(A_u,n_u,tau,order,packed)
  ##
  ## Computes the U* matrix at time tau in terms of A_u
  ## n_u : Number of system inputs
  ## If tau is a vector, computes U* at each tau and puts into a row vector:
  ## If packed==1
  ##     Ustar = [Ustar(tau_1) Ustar(tau_2) ...]
  ## else Ustar = [Ustar(tau_1); Ustar(tau_2) ...]
  ## Copyright (C) 1999 by Peter J. Gawthrop
  ## 	$Id$	

  if nargin<2
    n_u = 1;
  endif
  
  if nargin<3
    tau = 0;
  endif
  
  if nargin<4
    order = 0;
  endif
  
  if nargin<5
    packed=1;
  endif
  

  [n,m] = size(A_u);		# Size of composite A_u matrix
  N = m;			# Number of U* functions per input  
  nm = n/m;

  if (nm != n_u)&&(n!=m)	# Check consistency
    error("A_u must be square or be a column of square matrices");
  endif

  u_0 = ones(N,1);

  Ustar = [];
  for t = tau;
    ustar = [];
    for i = 1:n_u
      A_i = ppp_extract(A_u,i);
      Ak = A_i^order;
      eA = expm(A_i*t);
      if (packed==1)
	ustar = [ustar; zeros(1,(i-1)*N), (Ak*eA*u_0)', \
		 zeros(1,(n_u-i)*N)];
      else
	ustar = [ustar, (Ak*eA*u_0)'];
      endif
    endfor

    if (packed==1)
      Ustar = [Ustar ustar];
    else
      Ustar = [Ustar; ustar];
    endif
  endfor


endfunction