Overview
Comment:Put in View option to momitor progress graphically
Downloads: Tarball | ZIP archive
Timelines: family | ancestors | descendants | both | origin/master | trunk
Files: files | file ages | folders
SHA3-256: 323dea686653155aef3be942d82fabdb65d28f377a09c34d7ecf5c7323af469a
User & Date: gawthrop@users.sourceforge.net on 1999-12-07 23:09:56.000
Other Links: branch diff | manifest | tags
Context
1999-12-08
02:04:46
Removed bug - uj := MTTU(j,1); not commented out check-in: 77f5b6a386 user: gawthrop@users.sourceforge.net tags: origin/master, trunk
1999-12-07
23:09:56
Put in View option to momitor progress graphically check-in: 323dea6866 user: gawthrop@users.sourceforge.net tags: origin/master, trunk
1999-12-03
00:04:50
Version to eric.
-stdin switch added check-in: d32fbc4a33 user: gawthrop@users.sourceforge.net tags: origin/master, trunk
Changes
Unified Diff Ignore Whitespace Patch
Modified mttroot/mtt/bin/trans/m/mtt_optimise.m from [80100ee5cc] to [7de5a5a949]. 
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 = 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 = 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
MTT: Model Transformation Tools
[ GitHub | SourceHut | Sourceforge | Fossil RSS ]