#! /bin/sh
######################################
##### Model Transformation Tools #####
######################################
# Bourne shell script: dae2cse_r
# Differential-algebraic equations to constrained-state equations
# P.J.Gawthrop 14 June 1991, 8 Aug 1991, 2 April 1992, 14 April 1994, 28 Dec 94
# Copyright (c) P.J.Gawthrop 1991, 1992, 1994.
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.28 2005/09/07 17:04:12 geraint
## Fixes chi for cse representation.
##
## Revision 1.27 2005/03/21 11:50:39 gawthrop
## Don't write an empty cseo file
##
## Revision 1.26 2002/10/28 23:31:21 gawthrop
## Added additional term to MTTEdx to account for zdot terms on
## Right-Hand Side state equations
##
## Revision 1.25 2002/09/18 19:46:48 geraint
## Refined test to reduce shell loop overhead for unnecessary E matrix.
## Removed unnecessary attempt to optimise cse.
## Fixes sm representation.
##
## Revision 1.24 2002/09/16 08:08:00 geraint
## Merged changes from global-optimisation branch.
##
## Revision 1.23.2.3 2002/09/12 18:50:50 geraint
## Uncommented cse optimisations - they seem to work ok.
##
## Revision 1.23.2.2 2002/09/10 23:24:19 geraint
## Rationalised local and global optimisations.
## Fixes presentation of locally optimised code (ode view).
## Much more elegant :-)
##
## Revision 1.23.2.1 2002/09/03 23:44:43 geraint
## adding global optimisation (-optg).
##
## Revision 1.23 2002/07/10 17:43:05 geraint
## Added feature [ 562453 ] Optimisation of algebraic equations.
##
## Revision 1.22 2002/06/28 15:35:47 geraint
## Commented out aej.r generation (not used yet).
##
## Revision 1.21 2002/06/28 10:13:40 geraint
## Includes fix_c.r in ese2rdae and def2write_r to eliminate occurrances of x**y.
##
## Revision 1.20.2.1 2002/06/05 11:14:50 geraint
## ae.r now generated using def2write_r like cse?.r
## fix_c.r called at ese2rdae stage so that pow gets fixed in ae.r.
##
## Revision 1.20 2002/04/28 18:41:26 geraint
## Fixed [ 549658 ] awk should be gawk.
## Replaced calls to awk with call to gawk.
##
## Revision 1.19 2001/10/26 01:01:49 geraint
## fixcc when rdae_is_dae (-cr).
##
## Revision 1.18 2001/10/05 23:37:32 geraint
## Fixed assignment statement in ae.r when RHS=0.
##
## Revision 1.17 2001/07/27 23:29:10 geraint
## Optimises only when requested (-opt).
##
## Revision 1.16 2001/07/13 04:54:04 geraint
## Branch merge: numerical-algebraic-solution back to main.
##
## Revision 1.15.2.4 2001/06/26 00:55:48 geraint
## Writes algebraic equation Jacobian _aej.r (not used yet).
##
## Revision 1.15.2.3 2001/05/09 00:19:22 geraint
## Fixed EOF error when MTTNYZ=0.
##
## Revision 1.15.2.2 2001/05/05 20:50:16 geraint
## Fixed errors when MTTNx=0.
##
## Revision 1.15.2.1 2001/05/04 04:07:24 geraint
## Numerical solution of algebraic equations.
## sys_ae.cc written for unsolved inputs.
## Solution of equations using hybrd from MINPACK (as used by Octave fsolve).
##
## Revision 1.15 2001/03/19 02:28:52 geraint
## Branch merge: merging-ode2odes-exe back to MAIN.
##
## Revision 1.14.2.1 2001/03/19 00:29:08 geraint
## Parse switches (-A) before calling def2write_r.
## Update $1_def.* instead of removing.
##
## Revision 1.14 2000/12/28 12:24:35 peterg
## *** empty log message ***
##
## Revision 1.13 2000/10/11 08:52:46 peterg
## Creates csex (cse with dxe only) rep.
##
## Revision 1.12 2000/10/10 21:00:58 peterg
## New code genration
##
## Revision 1.11 1998/11/26 09:18:55 peterg
## Incluse subs.r
##
## Revision 1.10 1998/11/18 13:50:29 peterg
## Removed writeing of EYz matrix
##
## Revision 1.9 1998/11/18 10:53:38 peterg
## Put in some more "IF MTTNx>0 THEN" to avoid error messages when no
## states.
##
## Revision 1.8 1998/11/10 08:54:34 peterg
## Put in "IF MTTNx>0 THEN" to prevent probs when Nx=0
## -- still a couple of apparent error messages - but answers now
## correct
##
## Revision 1.7 1998/10/05 10:46:15 peterg
## Commented out redundant MTTY := MTTY + MTTEyx*MTTEdX;
##
## Revision 1.6 1998/07/19 12:44:35 peterg
## Set MTTYz := 0 if the array is empty - avoids irritating error
## message.
##
## Revision 1.5 1998/05/20 15:23:26 peterg
## Put MTTYz := MTTYz outsise the BEGIN/END
##
## Revision 1.4 1998/05/20 15:13:09 peterg
## Writes out algebraic equations (if any).
##
## Revision 1.3 1998/03/03 09:02:46 peterg
## Replaced MTTEyx*MTTEdX + MTTEyu*MTTdu; term
##
## Revision 1.2 1997/08/26 08:22:36 peterg
## Changed
## MTTY := MTTY + MTTEyx*MTTdX + MTTEyu*MTTdu;
## to
## MTTY := MTTY + MTTEyx*MTTEdX + MTTEyu*MTTdu;
##
## This sorts out the problem when dz appears in the output equation.
##
## Revision 1.1 1997/08/26 08:20:18 peterg
## Initial revision
##
## Revision 1.2 1996/08/25 09:57:30 peter
## Sorted out bug when MTTNz=0
##
## Revision 1.1 1996/08/15 16:47:02 peter
## Initial revision
##
###############################################################
#Explicit solution option and optimise option
solve=0; solve_msg=''
optimise=''; optimise_msg=''
while [ -n "`echo $1 | grep '^-'`" ]; do
case $1 in
-A )
solve=1
solve_msg=' with explicit solution of algebraic equations' ;;
-fixcc )
fixcc='-fixcc'
include=`echo 'in "'$MTT_LIB'/reduce/fix_c.r";'`
fix_msg='fixing c and cc code';
;;
-optimise_global )
optimise='-optimise_global'
optimise_msg=' with global optimisation' ;;
-optimise_local )
optimise='-optimise_local'
optimise_msg=' with local optimisation' ;;
*)
echo "$1 is an invalid argument - ignoring" ;;
esac
shift
done
# Create the reduce output code
def2write_r $fixcc $1 ae
def2write_r $fixcc $1 cse
def2write_r $fixcc $1 csex # Version without E matrix
def2write_r $fixcc $1 cseo
echo "Creating $1_ae.r $optimise_msg"
echo "Creating $1_cse.r $solve_msg $optimise_msg $fix_msg"
echo "Creating $1_csex.r $optimise_msg"
echo "Creating $1_cseo.r $optimise_msg"
# Remove the old log file
rm -f dae2cse_r.log
# Remove some files
rm -f $1_ae.r? $1_cse.r? $1_cseo.r?
# Use reduce to accomplish the transformation
$SYMBOLIC >dae2cse_r.log << EOF
%Read the formatting function
in "$MTTPATH/trans/reduce_matrix.r";
OFF Echo;
OFF Nat;
ON NERO;
%Fix c code if required
$include
in "$1_def.r";
MTTdxs := MTTdX; %Save the symbolic form of dX
in "$1_subs.r";
in "$1_dae.r";
%Create F_x, F_y matrices - assumming equations are
% linear in dZ
IF MTTNz>0 THEN
BEGIN
IF MTTNx>0 THEN
BEGIN
% Find MTTFx;
write "% Find MTTFx;";
matrix MTTFx(MTTNx,MTTNz);
FOR j := 1:MTTNz DO
BEGIN
dzj := MTTdZ(j,1);
FOR i := 1:MTTNx DO
MTTFx(i,j) := df(MTTdX(i,1), dzj, 1);
END;
END;
% Find MTTFy;
write "% Find MTTFy;";
matrix MTTFy(MTTNy,MTTNz);
FOR j := 1:MTTNz DO
BEGIN
dzj := MTTdZ(j,1);
FOR i := 1:MTTNy DO
MTTFy(i,j) := df(MTTy(i,1), dzj, 1);
END;
%Create G_x, G_u matrices
write "%Create G_x, G_u matrices ";
% Find MTTGx;
IF MTTNx>0 THEN
BEGIN
write "% Find MTTGx;";
matrix MTTGx(MTTNz,MTTNx);
FOR j := 1:MTTNx DO
BEGIN
xj := MTTX(j,1);
FOR i := 1:MTTNz DO
MTTGx(i,j) := df(MTTZ(i,1), xj, 1);
END;
END;
% Find MTTGu;
write "% Find MTTGu;";
matrix MTTGu(MTTNz,MTTNu);
FOR j := 1:MTTNu DO
BEGIN
uj := MTTu(j,1);
FOR i := 1:MTTNz DO
MTTGu(i,j) := df(MTTZ(i,1), uj, 1);
END;
%Create E matrices
write "%Create E matrices";
IF MTTNx>0 THEN
BEGIN
matrix MTTExx(MTTNx,MTTNx); MTTExx := MTTFx*MTTGx;
matrix MTTExu(MTTNx,MTTNu); MTTExu := MTTFx*MTTGu;
matrix MTTEyx(MTTNy,MTTNx); MTTEyx := MTTFy*MTTGx;
matrix MTTE(MTTNx,MTTNx); MTTE := MTTI - MTTExx;
END;
matrix MTTEyu(MTTNy,MTTNu); MTTEyu := MTTFy*MTTGu;
%% The following gets rid of the dZs; there must be a better way.
MTTdZ1 := 0;
MTTdZ2 := 0;
MTTdZ3 := 0;
MTTdZ4 := 0;
MTTdZ5 := 0;
MTTdZ6 := 0;
MTTdZ7 := 0;
MTTdZ8 := 0;
MTTdZ9 := 0;
MTTdZ10 := 0;
MTTdZ11 := 0;
MTTdZ12 := 0;
MTTdZ13 := 0;
MTTdZ14 := 0;
MTTdZ15 := 0;
MTTdZ16 := 0;
MTTdZ17 := 0;
MTTdZ18 := 0;
MTTdZ19 := 0;
IF MTTNx>0 THEN
BEGIN
MTTEdX := MTTdX; %Ie MTTEdX is MTTdX with the dz terms deleted ie EdX.
MTTdX := MTTdXs; %Restore the symbolic dX
%% Add on input derivative terms
MTTEdX := MTTEdX + MTTExu*MTTdu;
%% Add on output derivative terms
MTTEdx := MTTEdX + MTTEyx*(MTTE^(-1))*MTTEdX;
END;
%%%%%MTTY := MTTY + MTTEyx*MTTEdX;
%%% This causes the matrix mismatch
%%% MTTdXs and MTTdu need setting in _def.r file
MTTY := MTTY + MTTEyu*MTTdu;
IF MTTNx>0 THEN
MTTY := MTTY + MTTEyx*(MTTE^(-1))*MTTEdX;
END; %%of MTTNz>0
IF MTTNz=0 THEN
BEGIN
MTTEdX := MTTdX;
MTTE := MTTI;
END;
IF (MTTNyz>0) AND ($solve>0) THEN
BEGIN
%%%% Try and solve algebraic loops!!
%Create list of the relevant equations
MTT_eqns := {};
FOR i := 1:MTTNyz DO
MTT_eqns := append(MTT_eqns,{MTTyz(i,1)});
%Create list of the relevant unknowns
MTT_unknowns := {};
FOR i := 1:MTTNyz DO
MTT_unknowns := append(MTT_unknowns,{MTTUi(i,1)});
%Solve the algebraic equations symbolically
MTT_sol := solve(MTT_eqns,MTT_unknowns);
%The result seems to be in an extra list - I dont know why
% So remove the outer list with first.
% But only if more than one list element!
if MTTNyz>1 THEN
MTT_sol := first(MTT_sol);
%Substitute back into the equations
FOR i := 1:MTTNyz DO
BEGIN
MTT_sol_i := first(MTT_sol); MTT_sol := rest(MTT_sol);
set(lhs(MTT_sol_i),rhs(MTT_sol_i));
END;
% No algebraic variables left!
MTTNYz := 0;
END; % IF MTTNyz>0 and $solve
%OUT "$1_aej.r";
%IF (MTTNyz>0) THEN % as above
%BEGIN
% WRITE "MATRIX MTTyzj(",MTTNyz,",",MTTNyz,")";
% WRITE "%File: $1_aej.r";
% FOR i := 1:MTTNyz DO
% FOR j := 1:MTTNyz DO
% BEGIN
% didj := df(MTTyz(i,1),mkid('mttui,j));
% IF (didj NEQ 0) THEN
% WRITE "MTTyzj(",i,",",j,") := ",didj," +0";
% END;
%END;
%WRITE ";END;";
%SHUT "$1_aej.r";
IF MTTNyz>0 THEN % not $solve or solution failed
BEGIN
OUT "$1_ae.r1";
write "MATRIX MTTYZ(", MTTNyz, ",", 1, ")$";
SHUT "$1_ae.r1";
END;
OUT "$1_ae.r2";
write "%File: $1_ae.r";
in ("$1_ae_write.r");
write "END;";
SHUT "$1_ae.r2";
% Create the matrix declarations
OUT "$1_cse.r1";
write "%";
IF (MTTNx > 0) THEN
BEGIN
write "MATRIX MTTEdx(", MTTNx, ",", 1, ")$";
write "MATRIX MTTE(", MTTNx, ",", MTTNx, ")$";
END;
SHUT "$1_cse.r1";
OUT "$1_csex.r1";
write "%File:$1_csex.r1";
IF (MTTNx > 0) THEN
write "MATRIX MTTEdx(", MTTNx, ",", 1, ")$";
SHUT "$1_csex.r1";
IF MTTNy>0 THEN
BEGIN
OUT "$1_cseo.r1";
write "MATRIX MTTY(", MTTNy, ",", 1, ")$";
SHUT "$1_cseo.r1";
END;
%%Create the _cse.r file
OUT "$1_cse.r2";
write "%File: $1_cse.r";
in ("$1_cse_write.r");
write "in ""$1_cseo.r"";";
write "END;";
% % State derivative
% MTT_Matrix := MTTEdX$
% MTT_Matrix_name := "MTTEdX"$
% MTT_Matrix_n := MTTNx$
% MTT_Matrix_m := 1$
% Reduce_Matrix()$
% % Output
% MTT_Matrix := MTTY$
% MTT_Matrix_name := "MTTY"$
% MTT_Matrix_n := MTTNy$
% MTT_Matrix_m := 1$
% Reduce_Matrix()$
% % Inputs
% MTT_Matrix := MTTU$
% MTT_Matrix_name := "MTTU"$
% MTT_Matrix_n := MTTNu$
% MTT_Matrix_m := 1$
% Reduce_Matrix()$
% MTT_Matrix := MTTdU$
% MTT_Matrix_name := "MTTdU"$
% MTT_Matrix_n := MTTNu$
% MTT_Matrix_m := 1$
% Reduce_Matrix()$
% % E matrix
% MTT_Matrix := MTTE$
% MTT_Matrix_name := "MTTE"$
% MTT_Matrix_n := MTTNx$
% MTT_Matrix_m := MTTNx$
% Reduce_Matrix()$
% % Eyx matrix
% MTT_Matrix := MTTEyx$
% MTT_Matrix_name := "MTTEyx"$
% MTT_Matrix_n := MTTNy$
% MTT_Matrix_m := MTTNx$
% %Reduce_Matrix()$
% % Yz
% MTT_Matrix := MTTYz$
% MTT_Matrix_name := "MTTYz"$
% MTT_Matrix_n := MTTNyz$
% MTT_Matrix_m := 1$
% Reduce_Matrix()$
write ";END;";
SHUT "$1_cse.r2";
OUT "$1_csex.r2";
write "%File: $1_cse.r";
in ("$1_csex_write.r");
write "END;";
SHUT "$1_csex.r2";
%Write out the output equations
OUT "$1_cseo.r2";
write "%File: $1_cseo.r";
IF MTTNy>0 THEN
in ("$1_cseo_write.r");
write "END;";
SHUT "$1_cseo.r2";
quit;
EOF
touch $1_ae.r1 $1_ae.r2
touch $1_cseo.r1
touch $1_cseo.r2
cat $1_ae.r1 $1_ae.r2 > $1_ae.r
cat $1_cse.r1 $1_cse.r2 > $1_cse.r
cat $1_csex.r1 $1_csex.r2 > $1_csex.r
cat $1_cseo.r1 $1_cseo.r2 > $1_cseo.r
if [ ${optimise:-""} = "-optimise_global" ]; then
mtt_optimise global $1 ae
mtt_optimise global $1 cseo
mtt_optimise global $1 csex
elif [ ${optimise:-""} = "-optimise_local" ]; then
mtt_optimise local $1 ae
mtt_optimise local $1 cseo
mtt_optimise local $1 csex
fi
if [ "$solve" = "1" ]; then
echo "Setting MTTNyz=0 in $1_def.r and updating other $1_def files"
gawk '{
if ($1=="MTTNyz")
print "MTTNyz := 0;"
else print $0
}' $1_def.r > mtt_junk
# Make sure it preserves the time stamp!!
# and remove dependent reps
touch -r $1_def.r mtt_junk
mv mtt_junk $1_def.r
for file in `ls $1_def.*`; do
if [ $file != $1_def.r ]; then
ext=`echo $file|gawk -F\. '{print $2}'`
mtt -q $1 def $ext
fi
done
fi
# Now invoke the standard error handling.
mtt_error_r dae2cse_r.log