Overview
Comment:Addes Q weighting matrix
Downloads: Tarball | ZIP archive | SQL archive
Timelines: family | ancestors | descendants | both | origin/master | trunk
Files: files | file ages | folders
SHA3-256: 8e162b30dce1b77d90e609e7a6c6bab48e40aec1643d4fdfe3ea3e4f69c70e6d
User & Date: gawthrop@users.sourceforge.net on 2002-05-13 16:01:10
Other Links: branch diff | manifest | tags
Context
2002-05-13
16:11:09
Various changes check-in: 079339165e user: gawthrop@users.sourceforge.net tags: origin/master, trunk
16:01:10
Addes Q weighting matrix check-in: 8e162b30dc 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: b97e3c9beb user: gawthrop@users.sourceforge.net tags: origin/master, trunk
Changes

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
function [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,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)  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);
  
  U = par(free(:,1));
  U_all = Par(free(:,1),:);
endfunction



|

|



|












>
>
>
>






|







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, Q, extras)

  ## usage:  [U, U_all, Error, Y,its ] = ppp_nlin(system_name,x_0,par_0,sim,us,w,free, Q, extras)
  ##
  ## 
  
  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,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
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)


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







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







|







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





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







>
>
>
>







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







|







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(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
      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);
      if extras.visual
	figure(3);
	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];







|

|







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,Q,extras);
      if extras.visual
	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
function [par,Par,Error,Y,iterations,x] = \
      ppp_optimise(system_name,x_0,par_0,simpar,u,y_0,free,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

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

|














>














>
>
>







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





  n_th = length(i_s);
  err_old = inf;
  err_old_old = inf;
  err = 1e50;
  reduction = inf;
  predicted_reduction = 0;
  par = par_0;







|












>
>
>
>







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







|

|
|







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 + 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 + 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 ]