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