Overview
Comment:Updated to account for new nonlinear ppp
Downloads: Tarball | ZIP archive | SQL archive
Timelines: family | ancestors | descendants | both | origin/master | trunk
Files: files | file ages | folders
SHA3-256: 6632bcb60e10755fbfaeff588255efd5623c0152a4b085fd91c8e60d1e52f968
User & Date: gawthrop@users.sourceforge.net on 2001-05-26 15:46:38
Other Links: branch diff | manifest | tags
Context
2001-05-26
18:36:43
Further modifications. Now works on rcPPP
-- next jobs:
add identification to ppp_nlin_sim
create real-time ppp_nlin_run check-in: 3f62cd4a52 user: gawthrop@users.sourceforge.net tags: origin/master, trunk
15:46:38
Updated to account for new nonlinear ppp check-in: 6632bcb60e user: gawthrop@users.sourceforge.net tags: origin/master, trunk
2001-05-24
07:48:17
Include artwork in the cbg.fig file check-in: 39adc84a7e user: gawthrop@users.sourceforge.net tags: origin/master, trunk
Changes
Hide Diffs Side-by-Side Diffs Ignore Whitespace Patch

Modified mttroot/mtt/bin/mttrc.csh from [1bee2798af] to [689359fe62]. 

1
2

3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19









20
21
22
23
24
25
26

1

2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35


-
+

















+
+
+
+
+
+
+
+
+








#!/bin/csh
## Automatically generated from bashrc on Tue Apr 10 14:02:38 BST 2001 - DO NOT EDIT
## 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
︙
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

46
47
48
49
50
51
52

53
54
55
56
57
58
59








-








## 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
## If then else format
##
## 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 ***
##
︙
71
72
73
74
75
76
77
78

79
80
81
82
83
84
85
86
87
88
89
90
91
92

93
94
95
96
97
98
99
100
101
102
103
104

79
80
81
82
83
84
85

86
87
88
89
90
91
92
93
94
95
96
97
98
99

100
101
102



103
104
105
106
107
108
109








-
+













-
+


-
-
-








## 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
## 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 $1 by the mtt base path, eg /usr/share/mtt/latest
## 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

#if [ -z "$MTT_BASE" ]; then
#  echo mttrc requires one argument: eg mttrc /usr/share/mtt/latest
#else
  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'
︙
155
156
157
158
159
160
161
162


163
164
165
166
167
168
169
170
171
172





173
174
175
176
177




178
179
180
181

182
183
184
185
186
187
188
189
190






191
192
193
194
195
196



197
198

199
200
201
202
203
204
205

160
161
162
163
164
165
166

167
168
169
170
171
172
173





174
175
176
177
178
179




180
181
182
183
184
185
186

187
188
189
190






191
192
193
194
195
196
197
198
199



200
201
202
203

204
205
206
207
208
209
210









-
+
+





-
-
-
-
-
+
+
+
+
+

-
-
-
-
+
+
+
+



-
+



-
-
-
-
-
-
+
+
+
+
+
+



-
-
-
+
+
+

-
+






-

  # 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
  #setenv MKOCTFILE $MTT_LIB/octave/mkoctfile # This for no optimisation
    setenv MKOCTFILE mkoctfile

  # ode2odes.exe stuff

    # local system

    setenv PLAT "i686-pc-linux-gnu"
#    PREFIX "/usr/local"
    setenv PREFIX "/usr"
    setenv GCCVERS "2.95.2"
    setenv SRCOCTAVE "/cvs/octave"
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"
#    PLAT="mips-sgi-irix6.5"
#    PREFIX="/usr/people/bevangp/GNU"
#    GCCVERS="2.95.2"
#    SRCOCTAVE="${PREFIX}/../build/octave-2.1.33"

    # include paths

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

    # library paths

#    LOCTAVE "-L${PREFIX}/lib/octave -loctave -lcruft -loctinterp"
    setenv LOCTAVE "-L${PREFIX}/lib/octave -loctave -lcruft -loctinterp"
    setenv LKPATHSEA "-L${SRCOCTAVE}/kpathsea -lkpathsea"
    setenv LREADLINE " -L${SRCOCTAVE}/readline -lreadline"
    setenv LSYSTEM "-ldl -lm -lncurses"
    setenv LF2C "-L${PREFIX}/lib/gcc-lib/${PLAT}/${GCCVERS} -lg2c"
#    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

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

    # setenved variables
    # 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 " "
#fi

Modified mttroot/mtt/bin/trans/make_ssim from [dcfc761ddf] to [3ba8b84f9a].

︙
11
12
13
14
15
16
17
18







19
20
21
22
23
24
25
26
27
28


29
30
31
32
33
34
35
36
37

38
39


40
41
42
43
44
45
46

11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

27
28
29
30
31
32
33

34
35
36
37
38
39
40
41
42
43

44
45

46
47
48
49
50
51
52
53
54









+
+
+
+
+
+
+

-







-
+
+








-
+

-
+
+








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;
  ny = ny2/2;
  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+2*ny);
    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+2*ny);
	y = mtt_data(:,2:2:1+ny2);
      endif
      y_par = [y_par, mtt_data(:,3:2:2+2*ny)];
      y_par = [y_par, mtt_data(:,3:2:2+ny2)];
      x = [];
    endfor
  endif

endfunction
EOF

}
︙

Modified mttroot/mtt/bin/trans/mtt_header from [dbbc0973db] to [63911c75bc].

︙
8
9
10
11
12
13
14



15
16
17
18
19
20
21

8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24








+
+
+








# 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
##
︙
284
285
286
287
288
289
290
291

292
293
294
295
296
297
298

287
288
289
290
291
292
293

294
295
296
297
298
299
300
301








-
+








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

Modified mttroot/mtt/lib/comp/compound/Sensitivity/sCS/sCS_lbl.txt from [9950653fbb] to [da775c33bd]. 

1
2
3
4
5
6
7
8
9



10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

31
32
33
34
35
36
37
38

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

33
34
35
36
37
38
39
40
41










+
+
+




















-
+









%% 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		lin	effort,c;c_s
	c		slin	effort,c;c_s

% Component type SS
	[in]	 SS	external,external

% Component type Se
	e_0		SS	e_0;e_0s

Modified mttroot/mtt/lib/control/PPP/ppp_nlin.m from [0eacb03d2f] to [46741c2bd0]. 

1

2
3

4
5
6
7
8
9
10
11
12
13


14
15
16

17
18
19
20
21
22
23

24
25
26
27
28
29
30


1
2

3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

18
19
20
21
22
23
24

25
26
27
28
29
30
31
32

-
+

-
+










+
+


-
+






-
+








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

  ## usage:  U = ppp_nlin (system_name,x_0,u,t,par_0,free,w,weight)
  ## 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_us,n_tau] = size(us);
   [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(system_name,x_0,us,t,par_0,free,w,extras);
  [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

Modified mttroot/mtt/lib/control/PPP/ppp_nlin_sim.m from [93a2ec563f] to [7a76a326ac]. 

1

2
3
4
5
6
7
8


1
2
3
4
5
6
7
8

-
+








function [y,x,u,t,U,U_c,U_l] = ppp_nlin_sim (system_name,A_u,tau,t_ol,N_ol,w,extras)
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
︙
18
19
20
21
22
23
24
25

26

27

28
29
30


31
32
33
34

35
36
37
38
39
40






41
42






43
44
45
46


47
48
49
50
51
52
53
54
55
56

57
58
59
60
61
62
63
64
65

66
67
68
69
70
71
72



73
74
75
76
77
78
79
80
81
82
83
84
85
86
87

88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107

108
109
110
111
112
113
114

18
19
20
21
22
23
24

25
26
27

28
29
30

31
32
33
34
35
36
37






38
39
40
41
42
43


44
45
46
47
48
49
50
51


52
53
54
55
56
57
58
59
60
61
62

63
64
65
66
67
68
69
70
71

72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96

97
















98
99
100

101
102
103
104
105
106
107
108








-
+

+
-
+


-
+
+




+
-
-
-
-
-
-
+
+
+
+
+
+
-
-
+
+
+
+
+
+


-
-
+
+









-
+








-
+







+
+
+














-
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-



-
+








  endif
  
  

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

  ## System details 
  ## 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;", 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
  ## 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 = 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);
  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
  T_ol = t_ol(n_t)+dt
  

  ## -- 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)]; 
  ws = w*ones(1,n_tau);
#   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_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_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
		  
  
  ## 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);
  [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;
︙
123
124
125
126
127
128
129
130

131
132

133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152

153
154
155

156
157
158


159
160
161
162
163
164
165
166
167
168




169
170

171
172
173

174
175
176
177
178
179







180
181
182
183



184
185
186
187



188
189
190


191
192
193
194
195
196
197
198
199
200
201
202

117
118
119
120
121
122
123

124
125

126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145

146

147

148
149
150

151
152
153
154
155
156
157
158




159
160
161
162
163

164

165

166
167





168
169
170
171
172
173
174
175



176
177
178
179



180
181
182
183


184
185
186
187




188
189












-
+

-
+



















-
+
-

-
+


-
+
+






-
-
-
-
+
+
+
+

-
+
-

-
+

-
-
-
-
-
+
+
+
+
+
+
+

-
-
-
+
+
+

-
-
-
+
+
+

-
-
+
+


-
-
-
-


-
-
-
-

  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;
  U = expm(-A_u*T_ol)*U;

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

    ## Compute U
    ## Compute U by optimisation
    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);
      [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)
  
    ##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
    u_ol = U'*u_star_t(1:n_U,:);

    ## Simulate system 
    ## Simulate system over one ol interval
    ## (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);
    ##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(:,1:n_t)];
    x = [x x_ol(:,1:n_t)];
    u = [u u_ol(:,1:n_t)];
    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];
    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); 
    t = [t; t_ol+t_last*ones(n_t,1) ];
    t_last = t_last + T_ol; 
  endfor

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

Modified mttroot/mtt/lib/control/PPP/ppp_optimise.m from [a4cd6bb37d] to [b7facdcf98].

︙
24
25
26
27
28
29
30



31
32
33
34
35
36
37

24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40








+
+
+








  ######################################
  
  ###############################################################
  ## 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
  ##
︙
51
52
53
54
55
56
57



58
59
60
61
62
63
64

54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70








+
+
+








    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;
︙
88
89
90
91
92
93
94










95
96
97
98
99
100
101

94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117








+
+
+
+
+
+
+
+
+
+








  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
︙

Modified mttroot/mtt/lib/control/PPP/ppp_ustar.m from [7e105f13fe] to [0be0149c8e]. 

1

2
3

4
5
6
7

8

9
10
11
12
13
14
15
16
17
18
19
20
21
22




23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

42





43


44




45
46
47
48


1
2

3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

49
50
51
52
53
54
55
56

57
58
59
60
61
62
63
64

-
+

-
+




+

+














+
+
+
+



















+
-
+
+
+
+
+

+
+
-
+
+
+
+





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

  ## usage:  Us = ppp_ustar(A_u,n_u,tau)
  ## 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)];
	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];
      Ustar = [Ustar ustar];
    else
      Ustar = [Ustar; ustar];
    endif
  endfor


endfunction

MTT: Model Transformation Tools
[ GitHub | SourceHut | Sourceforge | Fossil RSS ]