Overview
Comment:error --> err to avoid name clash with built in function
Downloads: Tarball | ZIP 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.000
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
Unified Diff 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

  ######################################
  
  ###############################################################
  ## Version control history
  ###############################################################
  ## $Id$
  ## $Log$



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








>
>
>








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


  [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;
  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("  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)&&\
	(iterations<extras.max_iterations)

    iterations = iterations + 1; # Increment iteration counter

    [y,y_par,x] = eval(sim_command); # Simulate
    [N_data,N_y] = size(y);









|
|
|













|












|








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);
  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", 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(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


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

    ##Some diagnostics
    Error = [Error error];	# 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("  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)));








|





|






|










|
|










|
|


|





|








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
    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
      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 = 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
	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
    err_old_old = err_old;	# Save old error
    err_old = err;		# Save error

    ##Some diagnostics
    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", 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 ]