1
2
3
4
5
6
7
8
|
function [theta,Theta,Error,Y,iterations] = mtt_optimise (system_name,y_s,theta_0,method,free,weight,criterion,max_iterations,alpha)
## Usage: [theta,Theta,Error,Y,iterations] = mtt_optimise (system_name,y_s,theta_0,method,free,weight,criterion,max_iterations,alpha)
## 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
## method "time" or "freq"
|
|
|
1
2
3
4
5
6
7
8
|
function [theta,Theta,Error,Y,iterations] = mtt_optimise (system_name,y_s,theta_0,method,free,weight,criterion,max_iterations,alpha,View)
## Usage: [theta,Theta,Error,Y,iterations] = mtt_optimise (system_name,y_s,theta_0,method,free,weight,criterion,max_iterations,alpha)
## 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
## method "time" or "freq"
|
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
|
endif
if nargin <8
max_iterations = 25;
endif
if nargin<9
alpha = 1.0;
endif
if (!strcmp(method,"time"))&&(!strcmp(method,"freq"))
error("method must be either time or freq")
endif
[n_data,n_y] = size(y_s);
n_th = length(free);
error_old = 1e20;
error=1e10;
theta = theta_0;
Theta = [];
Error = [];
Y = [];
iterations = -1;
while (abs(error_old-error)>criterion)&&(iterations<max_iterations)
iterations = iterations + 1;
error_old_old = error_old;
error_old = error;
eval(sprintf("[t,y,y_theta] = mtt_s%s(system_name,theta,free);",method)); # Simulate system
error = error
## Diagnostics
Error = [Error error]; # Save error
Theta = [Theta theta]; # Save parameters
Y = [Y y]; # Save output
## Update the estimate if we are not done yet.
if (abs(error_old-error)>criterion)&&(iterations<max_iterations)
if error>error_old+criterion # Halve step size and try again
factor = 10;
disp(sprintf("step/%i",factor));
error = error_old; # Go back to previous error value
error_old = inf; # Don't let it think its converged
theta(free) = theta(free) + step; # Reverse step
step = step/factor; # new stepsize
else # Recompute step size
tol = 1e-5;
step = pinv(JJ,tol)*J; # Step size
#step = pinv(JJ)*J; # Step size (built in tol)
endif
theta(free) = theta(free) - step; # Increment parameters
endif
endwhile
endfunction
|
|
>
>
>
>
|
|
>
|
>
>
>
>
>
>
>
|
|
>
|
|
|
>
|
|
|
|
|
|
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
|
endif
if nargin <8
max_iterations = 25;
endif
if nargin<9
alpha = 0.1;
endif
if nargin<10
View = 0;
endif
if (!strcmp(method,"time"))&&(!strcmp(method,"freq"))
error("method must be either time or freq")
endif
[n_data,n_y] = size(y_s);
n_th = length(free);
error_old = inf;
error=1e50;
theta = theta_0;
Theta = [];
Error = [];
Y = [];
iterations = -1;
while (abs(error_old-error)>criterion)&&(iterations<max_iterations)
iterations = iterations + 1;
error_old_old = error_old;
error_old = error;
eval(sprintf("[t,y,y_theta] = \
mtt_s%s(system_name,theta,free);",method)); # Simulate system
if View
xlabel("");
title(sprintf("mtt_optimise: Weighted actual and estimated Interation %i", iterations));
plot(t,y.*weight,t,y_s.*weight);
endif
error = 0;
J = zeros(n_th,1);
JJ = zeros(n_th,n_th);
for i = 1:n_y
E = weight(:,i).*(y(:,i) - y_s(:,i)); # Weighted error
error = error + (E'*E); # Sum the error
Weight = weight(:,i)*ones(1,n_th); # Sensitivity weight
y_theta_w = Weight.*y_theta(:,i:n_y:n_y*n_th); # Weighted sensitivity
J = J + real(y_theta_w'*E); # Jacobian
JJ = JJ + real(y_theta_w'*y_theta_w); # Newton Euler approx Hessian
endfor
error
if error<(error_old+criterion) # Save some diagnostics
Error = [Error error]; # Save error
Theta = [Theta theta]; # Save parameters
Y = [Y y]; # Save output
endif
## Update the estimate if we are not done yet.
if (abs(error_old-error)>criterion)&&(iterations<max_iterations)
if error>(error_old+criterion) # Reduce step size and try again
factor = 10;
disp(sprintf("%2.2f*step",alpha));
error = error_old; # Go back to previous error value
error_old = inf; # Don't let it think its converged
theta(free) = theta(free) + step; # Reverse step
step = alpha*step; # new stepsize
else # Recompute step size
tol = 1e-5;
step = pinv(JJ,tol)*J; # Step size
#step = pinv(JJ)*J; # Step size (built in tol)
endif
theta(free) = theta(free) - step; # Increment parameters
endif
## theta
endwhile
endfunction
|