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
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
|
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
+
-
+
|
function [par,Par,Error,Y,iterations] = ppp_optimise(system_name,x_0,u,t,par_0,free,y_0,extras);
## Usage: [par,Par,Error,Y,iterations] = ppp_optimise(system_name,x_0,u,t,par_0,free,y_0,weight,extras);
## system_name String containg system name
## y_s actual system output
## theta_0 initial parameter estimate
## free Indices of the free parameters within theta_0
## weight Weighting function - same dimension as y_s
## extras.criterion convergence criterion
## extras.max_iterations limit to number of iterations
## extras.v Initial Levenberg-Marquardt parameter
function [par,Par,Error,Y,iterations] = \
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] = 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.1 2000/12/28 11:58:07 peterg
## Put under CVS
##
###############################################################
## 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_y,n_data] = size(y_0);
[n_data,n_y] = size(y_0);
n_th = length(i_s);
error_old = inf;
error_old_old = inf;
error = 1e50;
reduction = 1e50;
par = par_0;
|
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
|
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
|
-
+
-
-
-
-
-
+
+
+
+
+
-
|
while (abs(reduction)>extras.criterion)&&\
(abs(error)>extras.criterion)&&\
(iterations<extras.max_iterations)
iterations = iterations + 1; # Increment iteration counter
[y,y_par] = eval(sprintf("%s_sim(x_0,u,t,par,i_s)", system_name)); # Simulate
[y,y_par] = eval(sim_command); # Simulate
##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
E = y(i,:) - y_0(i,:); # Error
error = error + (E*E'); # Sum the error over outputs
y_par_i = y_par(i:n_y:n_y*n_th,:);
J = J + y_par_i*E'; # Jacobian
JJ = JJ + y_par_i*y_par_i'; # Newton Euler approx Hessian (with LM
E = y(:,i) - y_0(:,i); # Error in ith output
error = error + (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
# term
|
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
115
116
117
118
119
120
121
|
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
|
-
-
-
-
-
+
-
-
-
-
-
-
|
printf(" ratio: %g\n", r);
printf(" L-M param: %g\n", v);
printf(" parameters: ");
for i_th=1:n_th
printf("%g ", par(i_t(i_th)));
endfor
printf("\n");
Y = [Y; y]; # Save output
Y = [Y y]; # Save output
|