Overview
Comment:Addes Q weighting matrix
Downloads: Tarball | ZIP archive
Timelines: family | ancestors | descendants | both | origin/master | trunk
Files: files | file ages | folders
SHA3-256: 302ba43d9cb69ad58c1eca69599d2c68b95b917f20e83b8602bdfee7920a4093
User & Date: gawthrop@users.sourceforge.net on 2002-05-13 16:01:10.000
Other Links: branch diff | manifest | tags
Context
2002-05-13
16:11:09
Various changes check-in: 12e17795c6 user: gawthrop@users.sourceforge.net tags: origin/master, trunk
16:01:10
Addes Q weighting matrix check-in: 302ba43d9c user: gawthrop@users.sourceforge.net tags: origin/master, trunk
08:52:23
FIRST parameter of simpar now specifies first printed point -
simulation still starts at zero check-in: eb24171cc0 user: gawthrop@users.sourceforge.net tags: origin/master, trunk
Changes
Side-by-Side Diff Ignore Whitespace Patch
Modified mttroot/mtt/lib/control/PPP/ppp_nlin.m from [543adc29f3] to [d0ab21bae9]. 
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


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

-
+

-
+



-
+












+
+
+
+






-
+








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

  ## usage:  [U, U_all, Error, Y,its ] = ppp_nlin (system_name,x_0,par_0,sim,us,w,free,extras)
  ## usage:  [U, U_all, Error, Y,its ] = ppp_nlin(system_name,x_0,par_0,sim,us,w,free, Q, extras)
  ##
  ## 
  
  if nargin<8
  if nargin<9
    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);
  
  if nargin<8
    Q = ones(n_y,1);
  endif
 
  ## Checks
  if (n_us<>n_u)
    error(sprintf("Inputs (%i)  differenct to system inputs (%i)", n_us, n_u));
  endif
    
  ##Optimise
  [par,Par,Error,Y,its] = ppp_optimise(s_system_name,x_0,par_0,sim,us,w,free,extras);
  [par,Par,Error,Y,its] = ppp_optimise(s_system_name,x_0,par_0,sim,us,w,free,Q,extras);
  
  U = par(free(:,1));
  U_all = Par(free(:,1),:);
endfunction
Modified mttroot/mtt/lib/control/PPP/ppp_nlin_run.m from [a7eaf434e3] to [37eb564d85]. 
1

2
3
4
5
6
7
8


1
2
3
4
5
6
7
8

-
+








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)
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,Q,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
︙
27
28
29
30
31
32
33
34

35
36
37
38
39
40
41

27
28
29
30
31
32
33

34
35
36
37
38
39
40
41








-
+








  ## 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
  if nargin<8
    extras.alpha = 0.1;
    extras.criterion = 1e-5;
    extras.emulate_timing = 0;
    extras.max_iterations = 10;
    extras.simulate = 1;
    extras.v = 1e-5;
    extras.verbose = 0;
︙
57
58
59
60
61
62
63




64
65
66
67
68
69
70

57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74








+
+
+
+








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

  if nargin<8
    Q = ones(n_y,1);
  endif
 
  ## 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));
  x_0_models = x_0s;

  ## Times
︙
161
162
163
164
165
166
167
168

169
170
171
172
173
174
175

165
166
167
168
169
170
171

172
173
174
175
176
177
178
179








-
+








	
	## Optimise
	tick = time;
	[pars,Par,Error,Y,its] = \
	    ppp_optimise(s_system_name,x_0_models,pars,simpar_est,u_star_t,y_est,i_par,extras);
	
	if extras.visual
	  figure(2);
	  figure(11);
	  title("Parameter optimisation"); 
	  II = [1:length(y_est)]; plot(II,y_est,"*", II,Y);
	endif
	
	est_time = time-tick;  
	t_est = [t_est;est_time];
	its_est = [its_est; its-1];
︙
201
202
203
204
205
206
207
208

209
210

211
212
213
214
215
216
217

205
206
207
208
209
210
211

212
213

214
215
216
217
218
219
220
221








-
+

-
+








      x_next = x_next(n_t+1,:)'; # Initial state for next time
      x_nexts(1:2:(2*n_x)-1) = x_next; # And for internal sensitivity model
      
      ## 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);
      [U, U_all, Error, Y, its] = ppp_nlin(system_name,x_nexts,pars,simpars,u_star_tau,w_s,i_ppp,Q,extras);
      if extras.visual
	figure(3);
	figure(12);
	title("PPP optimisation");
	II = [1:length(w_s)]; plot(II,w_s,"*", II,Y);
	figure(1);
	endif

      ppp_time = time-tick;  
      t_ppp = [t_ppp;ppp_time];
︙
Modified mttroot/mtt/lib/control/PPP/ppp_optimise.m from [fc79616ba1] to [bf1568c82a]. 
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

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


-
+














+














+
+
+








function [par,Par,Error,Y,iterations,x] = \
      ppp_optimise(system_name,x_0,par_0,simpar,u,y_0,free,extras);
      ppp_optimise(system_name,x_0,par_0,simpar,u,y_0,free,Q,extras);
  ## Levenberg-Marquardt optimisation for PPP/MTT
  ## Usage: [par,Par,Error,Y,iterations,x] = ppp_optimise(system_name,x_0,par_0,simpar,u,y_0,free[,extras]);
  ##  system_name     String containing system name
  ##  x_0             Initial state
  ##  par_0           Initial parameter vector estimate
  ##  simpar          Simulation parameters:
  ##        .first      first time
  ##        .dt         time increment
  ##        .stepfactor Euler integration step factor
  ##  u               System input (column vector, each row is u')
  ##  y_0             Desired model output
  ##  free            one row for each adjustable parameter
  ##                  first column parameter indices
  ##                  second column corresponding sensitivity indices
  ##  Q               vector of positive output weights.
  ##  extras (opt)    optimisation parameters
  ##        .criterion convergence criterion
  ##        .max_iterations limit to number of iterations
  ##        .v  Initial Levenberg-Marquardt parameter

  ###################################### 
  ##### Model Transformation Tools #####
  ######################################
  
  ###############################################################
  ## Version control history
  ###############################################################
  ## $Id$
  ## $Log$
  ## Revision 1.9  2002/05/08 10:14:21  gawthrop
  ## Idetification now OK (Moved data range in ppp_optimise by one sample interval)
  ##
  ## Revision 1.8  2002/04/23 17:50:39  gawthrop
  ## error --> err to avoid name clash with built in function
  ##
  ## Revision 1.7  2001/08/10 16:19:06  gawthrop
  ## Tidied up the optimisation stuff
  ##
  ## Revision 1.6  2001/07/03 22:59:10  gawthrop
︙
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

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








-
+












+
+
+
+








  ## Copyright (C) 1999,2000 by Peter J. Gawthrop
  sim_command = sprintf("%s_ssim(x_0,par,simpar,u,i_s)", system_name);

  ## Extract indices
  i_t = free(:,1);		# Parameters
  i_s = free(:,2)';		# Sensitivities

  if nargin<8
  if nargin<9
    extras.criterion = 1e-5;
    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

  if nargin<8
    Q = ones(n_y,1);
  endif
  
  n_th = length(i_s);
  err_old = inf;
  err_old_old = inf;
  err = 1e50;
  reduction = inf;
  predicted_reduction = 0;
  par = par_0;
︙
126
127
128
129
130
131
132
133

134
135
136


137
138
139
140
141
142
143

134
135
136
137
138
139
140

141
142


143
144
145
146
147
148
149
150
151








-
+

-
-
+
+








    ##Evaluate error, cost derivative J and cost second derivative JJ
    err = 0; 
    J = zeros(n_th,1);
    JJ = zeros(n_th,n_th);
    
    for i = 1:n_y
      E = y(:,i) - y_0(:,i);	#  Error in ith output
      err = err + (E'*E);	# Sum the squared error over outputs
      err = err + Q(i)*(E'*E);	# Sum the squared error over outputs
      y_par_i = y_par(:,i:n_y:n_y*n_th); # sensitivity function (ith output)
      J  = J + y_par_i'*E;	# Jacobian
      JJ = JJ + y_par_i'*y_par_i; # Newton Euler approx Hessian
      J  = J + Q(i)*y_par_i'*E;	# Jacobian
      JJ = JJ + Q(i)*y_par_i'*y_par_i; # Newton Euler approx Hessian
    endfor

    if iterations>1 # Adjust the Levenberg-Marquardt parameter
      reduction = err_old-err;
      predicted_reduction =  2*J'*step + step'*JJ*step;
      r = predicted_reduction/reduction;
      if (r<0.25)||(reduction<0)
︙
MTT: Model Transformation Tools
[ GitHub | SourceHut | Sourceforge | Fossil RSS ]