# Acausal bond graph to causal bond graph: mfile format

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.12  2001/05/08 15:18:10  gawthrop
## Added trig and hyperbolic functions to argument exclusion list
##
## Revision 1.11  2000/09/19 11:14:30  peterg
## Now writes the first component type header correctely
##
## Revision 1.10  2000/01/26 10:11:10  peterg
## Added I component
##
## Revision 1.9  1999/11/10 00:47:08  peterg

  else {
    alias[++i] = $2;
    if (length(sensitivity)>0) {
      s_arg = sprintf("%ss",$2);
      alias[++i] = s_arg;
      s_arg = sprintf(";%s",s_arg);
      Component = substr(Component,2);
      lin = "slin"
    }
   else {
      lin = "lin"
    }
  
 }

    ## Table of components
    cr["SS"] = "SS";  arg["SS"] = "external,external";

    cr["Se"] = "SS";  arg["Se"] = sprintf("external%s", s_arg);
    cr["Sf"] = "SS";  arg["Sf"] = sprintf("external%s", s_arg);

    cr["De"] = "SS";  arg["De"] = "external";
    cr["Df"] = "SS";  arg["Df"] = "external";

    cr["R"] = "lin";  arg["R"]  = sprintf("flow,%s%s", $2, s_arg);
    cr["C"] = "lin";  arg["C"]  = sprintf("effort,%s%s", $2, s_arg);
    cr["I"] = "lin";  arg["I"]  = sprintf("flow,%s%s", $2, s_arg);
    cr["R"] = lin;  arg["R"]  = sprintf("flow,%s%s", $2, s_arg);
    cr["C"] = lin;  arg["C"]  = sprintf("effort,%s%s", $2, s_arg);
    cr["I"] = lin;  arg["I"]  = sprintf("flow,%s%s", $2, s_arg);

    cr["CS"] = "lin";  arg["CS"]  = sprintf("effort,%s;x_0%s", $2, s_arg);
    cr["IS"] = "lin";  arg["IS"]  = sprintf("flow,%s;x_0%s", $2, s_arg);
    cr["CS"] = "lin";  arg["CS"]  = sprintf("effort,%s;%s_x0%s", $2, $2, s_arg);
    cr["IS"] = "lin";  arg["IS"]  = sprintf("flow,%s;%s_x0%s", $2, $2, s_arg);
    cr["INTFS"] = "lin";  arg["INTFS"]  = sprintf("%s_x0%s", $2, s_arg);

    cr["TF"] = "lin";  arg["TF"]  = sprintf("flow,%s%s", $2, s_arg);
    cr["GY"] = "lin";  arg["GY"]  = sprintf("flow,%s%s", $2, s_arg);

    cr["AE"] = "lin";  arg["AE"]  = sprintf("%s%s", $2, s_arg);
    cr["AF"] = "lin";  arg["AF"]  = sprintf("%s%s", $2, s_arg);

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)
function [y,x,u,t,p,UU,UU_c,UU_l,t_ppp,t_est] = 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)
  ## usage: [y,x,u,t,p,UU,UU_c,UU_l,t_ppp,t_est] = ppp_nlin_sim (system_name,i_ppp,i_par,A_u,w_s,N_ol,extras) 
  ##
  ## 
  
  ## 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;
    extras.estimate = 1;
  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;

  [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 = [];

  p = [];

  t_last = 0;
  UU = [];
  UU_l =[];
  UU_c =[];
  
  t_ppp = [];
  t_est = [];

  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

    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;  
    ppp_time = time-tick;  
    printf("Optimisation %i took %i iterations and %2.2f sec\n", i, \
	   length(Error), opt_time);
    t_ppp = [t_ppp;ppp_time];
    
    ## Generate control
    u_ol = u_star_t*U;		# Not used - just for show

    ## Simulate system over one ol interval
    par(i_ppp(:,3)) = pars(i_ppp(:,1)); # Update the simulation ppp weights
    [y_ol,ys_ol,x_ol] = eval(sprintf("%s_ssim(x_0s, pars, simpar, u_star_t);", s_system_name));
    [y_ol,x_ol] = eval(sprintf("%s_sim(x_0, par, simpar, u_star_t);", system_name));


    ## Tune parameters/states
    if (extras.estimate==1)
      tick = time;
      par_est = pars(i_par(:,1));
      p = [p; par_est'];
      pars = ppp_optimise(s_system_name,x_0s,pars,simpar,u_star_t,y_ol,i_par,extras);
      est_time = time-tick;  
      t_est = [t_est;est_time];
    endif

    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

# -*-makefile-*-
# Makefile for representation nppp
# File nppp_rep.make

#Copyright (C) 2000 by Peter J. Gawthrop
#Copyright (C) 2000,2001 by Peter J. Gawthrop

all: $(SYS)_nppp.$(LANG)

$(SYS)_nppp.view: $(SYS)_nppp.ps
	echo Viewing $(SYS)_nppp.ps; ghostview $(SYS)_nppp.ps&

$(SYS)_nppp.ps: $(SYS)_ode2odes.out s$(SYS)_ode2odes.out \