Overview
Comment:Idetification now OK (Moved data range in ppp_optimise by one sample interval)
Downloads: Tarball | ZIP archive | SQL archive
Timelines: family | ancestors | descendants | both | origin/master | trunk
Files: files | file ages | folders
SHA3-256: 2584b3af97baf8e3590eed0286e358c7577f05dcad01201cb1d550452e677b82
User & Date: gawthrop@users.sourceforge.net on 2002-05-08 10:14:21
Other Links: branch diff | manifest | tags
Context
2002-05-08
11:37:30
Preliminary version of new style DIY rep for nppp check-in: fe8e70d7e1 user: gawthrop@users.sourceforge.net tags: origin/master, trunk
10:14:21
Idetification now OK (Moved data range in ppp_optimise by one sample interval) check-in: 2584b3af97 user: gawthrop@users.sourceforge.net tags: origin/master, trunk
2002-05-07
23:50:34
Preliminary support for Matlab dynamically linked shared objects:
invoke with: mtt -cc sys rep mexglx
ode2odes support is not yet included.
check-in: db3acb3909 user: geraint@users.sourceforge.net tags: origin/master, trunk
Changes

Modified mttroot/mtt/lib/control/PPP/ppp_RT_sim_compute.m from [1a61cd1008] to [d2604d4e33].

15
16
17
18
19
20
21
22
23
24
25
26

  ## Set up u_star
  u_star = ppp_ustar(A_u_sim,1,t,0,0,n_u-n_U);

  ## Simulate
  par(i_ppp_sim(:,3)) = U;		# Update the simulation ppp weights
  [y_sim,x] = eval(sprintf("%s_sim(x_0_sim, par, simpar_sim, u_star);", \
			     system_name_sim));
  x_0_sim  = x(n_t,:)';     # Extract state for next time
  u_sim = (u_star*U);
endfunction








|
|



15
16
17
18
19
20
21
22
23
24
25
26

  ## Set up u_star
  u_star = ppp_ustar(A_u_sim,1,t,0,0,n_u-n_U);

  ## Simulate
  par(i_ppp_sim(:,3)) = U;		# Update the simulation ppp weights
  [y_sim,x] = eval(sprintf("%s_sim(x_0_sim, par, simpar_sim, u_star);", \
			   system_name_sim));
  x_0_sim  = x(n_t,:)';		# Extract state for next time
  u_sim = (u_star*U);
endfunction

Modified mttroot/mtt/lib/control/PPP/ppp_nlin_run.m from [3008238ce6] to [a7eaf434e3].

35
36
37
38
39
40
41

42
43
44
45
46
47
48
49
50
51
52
    extras.alpha = 0.1;
    extras.criterion = 1e-5;
    extras.emulate_timing = 0;
    extras.max_iterations = 10;
    extras.simulate = 1;
    extras.v = 1e-5;
    extras.verbose = 0;

  endif

  ##Estimate if we have some adjustable parameters
  estimating_parameters = (length(i_par)>0)
  
  ## Names
  s_system_name = sprintf("s%s", system_name);

  ## System details -- defines simulation within ol interval
  par = eval(sprintf("%s_numpar;", system_name));
  simpar = eval(sprintf("%s_simpar;", system_name));







>



|







35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
    extras.alpha = 0.1;
    extras.criterion = 1e-5;
    extras.emulate_timing = 0;
    extras.max_iterations = 10;
    extras.simulate = 1;
    extras.v = 1e-5;
    extras.verbose = 0;
    extras.visual = 0;
  endif

  ##Estimate if we have some adjustable parameters
  estimating_parameters = (length(i_par)>0);
  
  ## Names
  s_system_name = sprintf("s%s", system_name);

  ## System details -- defines simulation within ol interval
  par = eval(sprintf("%s_numpar;", system_name));
  simpar = eval(sprintf("%s_simpar;", system_name));
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
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201





202
203
204
205
206
207
208

      ## Generate input to actual system
      u_star_t = ppp_ustar(A_u,1,t_ol',0,0,n_u-n_U);

      ## Tune parameters/states
      if (estimating_parameters==1)
	## Save up the estimated parameters
	par_est = pars(i_par(:,1))
	p = [p; par_est'];

	## Set up according to interval length
	if (T_ol>T_ol_0) ## Truncate data
	  simpar_est.last = T_ol_0;
	  y_est = y_ol(1:n_t+1,:);
	else
	  simpar_est.last = T_ol;
	  y_est = y_ol;
	endif

	simpar_pred.last = T_ol_0; # Predicted length of next interval
	pars(i_ppp(:,1)) = U_old; # Update the simulation ppp weights
	
	## Optimise
	tick = time;
	[pars,Par,Error,Y,its] = \
	    ppp_optimise(s_system_name,x_0_models,pars,simpar_est,u_star_t,y_est,i_par,extras);




II = [1:length(y_est)]; plot(II,y_est,"*", II,Y)


	est_time = time-tick;  
	t_est = [t_est;est_time];
	its_est = [its_est; its-1];
      endif

      ## Update internal model
      par(i_ppp(:,3)) = U_old; # Update the internal model ppp weights

      if (estimating_parameters==1)
	par(i_par(:,3)) = pars(i_par(:,1)); # Update the internal model params
      endif

      simpar_model.last = T_ol;
      [y_model,x_model] = eval(sprintf("%s_sim(x_0_model, par, simpar_model, \
 					       u_star_t);",system_name));

      x_0 = x_model(n_ol+1,:)';	# Initial state of next interval

      x_0_model = x_0;
      x_0_models(1:2:(2*n_x)-1) = x_0_model;

      ## Compute U by optimisation
      tick = time;

      ## Predict state at start of next interval
      par(i_ppp(:,3)) = U;
      [y_next,x_next] = eval(sprintf("%s_sim(x_0, par, simpar, \
					     u_star_t);",system_name));
      x_next = x_next(n_t+1,:)'; # Initial state for next time
      x_nexts(1:2:(2*n_x)-1) = x_next; # And for internal sensitivity model
      
      ## Optimize for next interval      
      U_old = U;		# Save previous value
      U = expm(A_u*T_ol)*U;	# Initialise from continuation trajectory
      pars(i_ppp(:,1)) = U;	# Put initial value of U into the parameter vector
      [U, U_all, Error, Y, its] = ppp_nlin(system_name,x_nexts,pars,simpars,u_star_tau,w_s,i_ppp,extras);







      ppp_time = time-tick;  
      t_ppp = [t_ppp;ppp_time];
      its_ppp = [its_ppp; its-1];

      ## Total execution time
      T_total = time - t_start;







|


















>
>
>
>
|
>
>

















>


















|
>
>
>
>
>







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
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221

      ## Generate input to actual system
      u_star_t = ppp_ustar(A_u,1,t_ol',0,0,n_u-n_U);

      ## Tune parameters/states
      if (estimating_parameters==1)
	## Save up the estimated parameters
	par_est = pars(i_par(:,1));
	p = [p; par_est'];

	## Set up according to interval length
	if (T_ol>T_ol_0) ## Truncate data
	  simpar_est.last = T_ol_0;
	  y_est = y_ol(1:n_t+1,:);
	else
	  simpar_est.last = T_ol;
	  y_est = y_ol;
	endif

	simpar_pred.last = T_ol_0; # Predicted length of next interval
	pars(i_ppp(:,1)) = U_old; # Update the simulation ppp weights
	
	## Optimise
	tick = time;
	[pars,Par,Error,Y,its] = \
	    ppp_optimise(s_system_name,x_0_models,pars,simpar_est,u_star_t,y_est,i_par,extras);
	
	if extras.visual
	  figure(2);
	  title("Parameter optimisation"); 
	  II = [1:length(y_est)]; plot(II,y_est,"*", II,Y);
	endif
	
	est_time = time-tick;  
	t_est = [t_est;est_time];
	its_est = [its_est; its-1];
      endif

      ## Update internal model
      par(i_ppp(:,3)) = U_old; # Update the internal model ppp weights

      if (estimating_parameters==1)
	par(i_par(:,3)) = pars(i_par(:,1)); # Update the internal model params
      endif

      simpar_model.last = T_ol;
      [y_model,x_model] = eval(sprintf("%s_sim(x_0_model, par, simpar_model, \
 					       u_star_t);",system_name));

      x_0 = x_model(n_ol+1,:)';	# Initial state of next interval
##      x_0 = x_model(n_ol-1,:)';	# Initial state of next interval
      x_0_model = x_0;
      x_0_models(1:2:(2*n_x)-1) = x_0_model;

      ## Compute U by optimisation
      tick = time;

      ## Predict state at start of next interval
      par(i_ppp(:,3)) = U;
      [y_next,x_next] = eval(sprintf("%s_sim(x_0, par, simpar, \
					     u_star_t);",system_name));
      x_next = x_next(n_t+1,:)'; # Initial state for next time
      x_nexts(1:2:(2*n_x)-1) = x_next; # And for internal sensitivity model
      
      ## Optimize for next interval      
      U_old = U;		# Save previous value
      U = expm(A_u*T_ol)*U;	# Initialise from continuation trajectory
      pars(i_ppp(:,1)) = U;	# Put initial value of U into the parameter vector
      [U, U_all, Error, Y, its] = ppp_nlin(system_name,x_nexts,pars,simpars,u_star_tau,w_s,i_ppp,extras);
      if extras.visual
	figure(3);
	title("PPP optimisation");
	II = [1:length(w_s)]; plot(II,w_s,"*", II,Y);
	figure(1);
	endif

      ppp_time = time-tick;  
      t_ppp = [t_ppp;ppp_time];
      its_ppp = [its_ppp; its-1];

      ## Total execution time
      T_total = time - t_start;

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

24
25
26
27
28
29
30



31
32
33
34
35
36
37
  ######################################
  
  ###############################################################
  ## 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







>
>
>







24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
  ######################################
  
  ###############################################################
  ## Version control history
  ###############################################################
  ## $Id$
  ## $Log$
  ## 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
  ## 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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
  ## 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<8







<







53
54
55
56
57
58
59

60
61
62
63
64
65
66
  ## 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<8
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135

    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
    
    ## Use the last part of the simulation to compare with data
    y = y(1+N_data-n_data:N_data,:);
    y_par = y_par(1+N_data-n_data:N_data,:);

    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







|
|

<
<
<
<
<







116
117
118
119
120
121
122
123
124
125





126
127
128
129
130
131
132

    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
    
    ## Use the last part of the simulation to compare with data
    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


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