Overview
Comment:error --> err to avoid name clash with built in function
Downloads: Tarball | ZIP archive | SQL archive
Timelines: family | ancestors | descendants | both | origin/master | trunk
Files: files | file ages | folders
SHA3-256: 1288c7a627cbcd4806976ef1ab1ad68a8d47425932a90caacab5561c1e83783b
User & Date: gawthrop@users.sourceforge.net on 2002-04-23 17:50:39
Other Links: branch diff | manifest | tags
Context
2002-04-25
09:56:27
\% --> % to avoid awk error message check-in: 833d750d78 user: gawthrop@users.sourceforge.net tags: origin/master, trunk
2002-04-23
17:50:39
error --> err to avoid name clash with built in function check-in: 1288c7a627 user: gawthrop@users.sourceforge.net tags: origin/master, trunk
17:46:05
_sim.m now returns time as third argument check-in: 864f048dbc user: gawthrop@users.sourceforge.net tags: origin/master, trunk
Changes
Hide Diffs Side-by-Side Diffs Ignore Whitespace Patch

Modified mttroot/mtt/lib/control/PPP/ppp_optimise.m from [87f475bb1d] to [293204b700].

︙
24
25
26
27
28
29
30



31
32
33
34
35
36
37

24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40








+
+
+








  ######################################
  
  ###############################################################
  ## Version control history
  ###############################################################
  ## $Id$
  ## $Log$
  ## 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
  ## Fixed problems with argument passing for CRs
  ##
  ## Revision 1.5  2001/06/06 07:54:38  gawthrop
  ## Further fixes to make nonlinear PPP work ...
  ##
  ## Revision 1.4  2001/05/26 15:46:38  gawthrop
︙
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

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
113
114








-
-
-
+
+
+













-
+












-
+









  [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);
  error_old = inf;
  error_old_old = inf;
  error = 1e50;
  err_old = inf;
  err_old_old = inf;
  err = 1e50;
  reduction = inf;
  predicted_reduction = 0;
  par = par_0;
  Par = par_0;
  step = ones(n_th,1);
  Error = [];
  Y = [];
  iterations = 0;
  v = extras.v;			# Levenverg-Marquardt parameter.
  r = 1;			# Step ratio

  if extras.verbose		# Diagnostics
    printf("Iteration: %i\n", iterations);
    printf("  error:  %g\n", error);
    printf("  error:  %g\n", err);
    printf("  reduction:  %g\n", reduction);
    printf("  prediction: %g\n", predicted_reduction);
    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");
  endif
  
  while (abs(reduction)>extras.criterion)&&\
	(abs(error)>extras.criterion)&&\
	(abs(err)>extras.criterion)&&\
	(iterations<extras.max_iterations)

    iterations = iterations + 1; # Increment iteration counter

    [y,y_par,x] = eval(sim_command); # Simulate
    [N_data,N_y] = size(y);
︙
120
121
122
123
124
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
154
155
156
157
158
159
160
161
162
163
164


165
166
167

168
169
170
171
172
173

174
175
176
177
178
179
180

123
124
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


154
155
156
157
158
159
160
161
162
163
164
165


166
167
168
169

170
171
172
173
174
175

176
177
178
179
180
181
182
183








-
+





-
+






-
+










-
-
+
+










-
-
+
+


-
+





-
+









    if extras.verbose		# Diagnostics
##      printf("y and y_0\n");
##      [y,y_0]
    endif
    
    ##Evaluate error, cost derivative J and cost second derivative JJ
    error = 0; 
    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
      error = error + (E'*E);	# Sum the squared error over outputs
      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 = error_old-error;
      reduction = err_old-err;
      predicted_reduction =  2*J'*step + step'*JJ*step;
      r = predicted_reduction/reduction;
      if (r<0.25)||(reduction<0)
	v = 4*v;
      elseif r>0.75
	v = v/2;
      endif

      if reduction<0		# Its getting worse
	par(i_t) = par(i_t) + step; # rewind parameter
	error = error_old;	# rewind error
	error_old = error_old_old; # rewind old error
	err = err_old;	# rewind error
	err_old = err_old_old; # rewind old error
	if extras.verbose
	  printf(" Rewinding ....\n");
	endif
      endif
    endif

    ## Compute step using pseudo inverse
    JJL = JJ + v*eye(n_th);	# Levenberg-Marquardt term
    step =  pinv(JJL)*J;	# Step size
    par(i_t) = par(i_t) - step; # Increment parameters
    error_old_old = error_old;	# Save old error
    error_old = error;		# Save error
    err_old_old = err_old;	# Save old error
    err_old = err;		# Save error

    ##Some diagnostics
    Error = [Error error];	# Save error
    Error = [Error err];	# Save error
    Par = [Par par];		# Save parameters
    Y = [Y y];			# Save output

    if extras.verbose		# Diagnostics
      printf("Iteration: %i\n", iterations);
      printf("  error:  %g\n", error);
      printf("  error:  %g\n", err);
      printf("  reduction:  %g\n", reduction);
      printf("  prediction: %g\n", predicted_reduction);
      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)));
︙

MTT: Model Transformation Tools
[ GitHub | SourceHut | Sourceforge | Fossil RSS ]