Overview
Comment:Restuctured to be more logical.
Data is now in columns to be compatible with MTT.
Downloads: Tarball | ZIP archive | SQL archive
Timelines: family | ancestors | descendants | both | origin/master | trunk
Files: files | file ages | folders
SHA3-256: 2a56dd249c1a7138405387d1b6aebeb17fa907f89dfe0b8f9a73427eb8ac23b5
User & Date: gawthrop@users.sourceforge.net on 2001-04-04 08:36:25
Other Links: branch diff | manifest | tags
Context
2001-04-04
10:05:38
Reresentation for system identification for ppp check-in: 7483f97d8c user: gawthrop@users.sourceforge.net tags: origin/master, trunk
08:36:25
Restuctured to be more logical.
Data is now in columns to be compatible with MTT. check-in: 2a56dd249c user: gawthrop@users.sourceforge.net tags: origin/master, trunk
01:03:01
Fixed elimination of 0 in xxx(0...) when ... does not start with a digit. check-in: 66c73b4221 user: geraint@users.sourceforge.net tags: origin/master, trunk
Changes
Hide Diffs Unified Diffs Ignore Whitespace Patch

Modified mttroot/mtt/lib/control/PPP/ppp_optimise.m from [aff901223d] to [55529a32e0]. 

1


2
3
4
5
6
7








8
9
10
11















12
13

14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

function [par,Par,Error,Y,iterations] = ppp_optimise(system_name,x_0,u,t,par_0,free,y_0,extras);


  ## Usage: [par,Par,Error,Y,iterations] = ppp_optimise(system_name,x_0,u,t,par_0,free,y_0,weight,extras);
  ##  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








  ##  extras.criterion convergence criterion
  ##  extras.max_iterations limit to number of iterations
  ##  extras.v  Initial Levenberg-Marquardt parameter
















  ## Copyright (C) 1999,2000 by Peter J. Gawthrop


  ## Extract indices
  i_t = free(:,1); # Parameters
  i_s = free(:,2)'; # Sensitivities

  if nargin<8
    extras.criterion = 1e-5;
    extras.max_iterations = 10;
    extras.v = 1e-5;
    extras.verbose = 0;
  endif
  

  [n_y,n_data] = size(y_0);

  n_th = length(i_s);
  error_old = inf;
  error_old_old = inf;
  error = 1e50;
  reduction = 1e50;
  par = par_0;

|
>
>
|
|
|
|
|
|
>
>
>
>
>
>
>
>
|
|
|

>
>
>
>
>
>
>
>
>
>
>
>
>
>
>


>




|







|








1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
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

function [par,Par,Error,Y,iterations] = \
      ppp_optimise(system_name,x_0,par_0,simpar,u,y_0,free,extras);
  ## Levenberg-Marquardt optimisation for PPP/MTT
  ## Usage: [par,Par,Error,Y,iterations] = ppp_optimise(system_name,x_0,par_0,simpar,u,y_0,free[,extras]);
  ##  system_name     String containing system name
  ##  x_0             Initial state
  ##  par_0           Initial parameter vector estimate
  ##  simpar          Simulation parameters:
  ##        .first      first time
  ##        .dt         time increment
  ##        .stepfactor Euler integration step factor
  ##  u               System input (column vector, each row is u')
  ##  y_0             Desired model output
  ##  free            one row for each adjustable parameter
  ##                  first column parameter indices
  ##                  second column corresponding sensitivity indices
  ##  extras (opt)    optimisation parameters
  ##        .criterion convergence criterion
  ##        .max_iterations limit to number of iterations
  ##        .v  Initial Levenberg-Marquardt parameter

  ###################################### 
  ##### Model Transformation Tools #####
  ######################################
  
  ###############################################################
  ## Version control history
  ###############################################################
  ## $Id$
  ## $Log$
  ## Revision 1.1  2000/12/28 11:58:07  peterg
  ## Put under CVS
  ##
  ###############################################################


  ## Copyright (C) 1999,2000 by Peter J. Gawthrop

  sim_command = sprintf("%s_ssim(x_0,par,simpar,u,i_s)", system_name)
  ## Extract indices
  i_t = free(:,1); # Parameters
  i_s = free(:,2)'; # Sensitivities

  if nargin<9
    extras.criterion = 1e-5;
    extras.max_iterations = 10;
    extras.v = 1e-5;
    extras.verbose = 0;
  endif
  

  [n_data,n_y] = size(y_0);

  n_th = length(i_s);
  error_old = inf;
  error_old_old = inf;
  error = 1e50;
  reduction = 1e50;
  par = par_0;
︙︙
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


  while (abs(reduction)>extras.criterion)&&\
	(abs(error)>extras.criterion)&&\
	(iterations<extras.max_iterations)

    iterations = iterations + 1; # Increment iteration counter

    [y,y_par] = eval(sprintf("%s_sim(x_0,u,t,par,i_s)", system_name)); # Simulate

    ##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
      error = error + (E*E');	# Sum the error over outputs
      y_par_i = y_par(i:n_y:n_y*n_th,:);
      J  = J + y_par_i*E';	# Jacobian
      JJ = JJ + y_par_i*y_par_i'; # Newton Euler approx Hessian (with LM
				# term
    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)








|







|
|
|
|
|
<








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


  while (abs(reduction)>extras.criterion)&&\
	(abs(error)>extras.criterion)&&\
	(iterations<extras.max_iterations)

    iterations = iterations + 1; # Increment iteration counter

    [y,y_par] = eval(sim_command); # Simulate

    ##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)
︙︙
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
115
116
117
118
119
120
121

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

  endwhile

endfunction













<









<

<





<







|
<



<
<
<
<
<

104
105
106
107
108
109
110

111
112
113
114
115
116
117
118
119

120

121
122
123
124
125

126
127
128
129
130
131
132
133

134
135
136






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

  endwhile

endfunction

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