25
26
27
28
29
30
31
32
33
34
35
36
37
38
|
######################################
###############################################################
## 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
|
>
>
>
|
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
|
######################################
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.10 2002/05/13 16:01:09 gawthrop
## Addes Q weighting matrix
##
## 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
|
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
|
[y,y_par,x] = eval(sim_command); # Simulate
[N_data,N_y] = size(y);
if (N_y!=n_y)
mess = sprintf("n_y (%i) in data not same as n_y (%i) in model", n_y,N_y);
error(mess);
endif
|
|
>
>
>
>
>
|
|
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
|
[y,y_par,x] = eval(sim_command); # Simulate
[N_data,N_y] = size(y);
if (N_y!=n_y)
mess = sprintf("n_y (%i) in data not same as n_y (%i) in model", n_y,N_y);
error(mess);
endif
if ( (N_data-n_data)<1 )
error(sprintf("y_0 (%i) must be shorter than y (%i)", n_data, N_data));
endif
## Use the last part of the simulation to compare with data
## And shift back by one data point
y = y(N_data-n_data:N_data-1,:);
y_par = y_par(N_data-n_data:N_data-1,:);
##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
|