Overview
Comment:added -ae option to select algebraic equation solution method.
Downloads: Tarball | ZIP archive | SQL archive
Timelines: family | ancestors | descendants | both | origin/numerical-algebraic-solution | trunk
Files: files | file ages | folders
SHA3-256: 97ad53b7470e78ba1a03e4c247d51a64a419c58a9dc42f3057e89c6cfd8c6bd6
User & Date: geraint@users.sourceforge.net on 2001-06-05 03:20:40
Other Links: branch diff | manifest | tags
Context
2001-06-05
03:37:15
*** empty log message *** check-in: 4c82a03e36 user: geraint@users.sourceforge.net tags: origin/numerical-algebraic-solution, trunk
03:20:40
added -ae option to select algebraic equation solution method. check-in: 97ad53b747 user: geraint@users.sourceforge.net tags: origin/numerical-algebraic-solution, trunk
2001-05-09
00:19:22
Fixed EOF error when MTTNYZ=0. check-in: 414839b0b6 user: geraint@users.sourceforge.net tags: origin/numerical-algebraic-solution, trunk
Changes

Modified mttroot/mtt/bin/mtt from [6ea445f6e4] to [af07966fc8].

12
13
14
15
16
17
18





19
20
21
22
23
24
25
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30







+
+
+
+
+







# Copyright (C) 2001 by Peter J. Gawthrop

###############################################################
## Version control history
###############################################################
## $Header$
## $Log$
## Revision 1.309.2.1  2001/05/04 04:07:24  geraint
## Numerical solution of algebraic equations.
## sys_ae.cc written for unsolved inputs.
## Solution of equations using hybrd from MINPACK (as used by Octave fsolve).
##
## Revision 1.309  2001/04/28 03:15:03  geraint
## Fixed comment (interfered with "mtt help representations").
##
## Revision 1.308  2001/04/25 22:17:45  geraint
## Fixed icd.txt2 dependency.
##
## Revision 1.307  2001/04/15 21:15:41  geraint
1086
1087
1088
1089
1090
1091
1092



1093
1094
1095
1096
1097
1098
1099
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107







+
+
+








# Default not verbose
verbose=' -s'

# Default integration method
integration_method=implicit;

# Default algebraic equation solver
algebraic_solver=Reduce_Solver;

# Default no info
info_switch=''

# Default use m, not oct files
m='m';

# Default use ps files
1150
1151
1152
1153
1154
1155
1156























1157
1158
1159
1160
1161
1162
1163
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194







+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+







			integration_method=rk4;
			mtt_switches="$mtt_switches rk4";
			;;
		    *)
			echo $1 is an unknown integration method - use euler, rk4 or implicit;
                        exit;;
		esac;;
	-ae )
                mtt_switches="$mtt_switches $1";
		case $2 in
		    fsolve | hybrd)
			mtt_switches="$mtt_switches $2";
			algebraic_solver=Hybrd_Solver;
			shift;
			;;
		    hooke)
			mtt_switches="$mtt_switches $2";
			algebraic_solver=HJ_Solver;
			shift;
			;;
		    reduce)
			mtt_switches="$mtt_switches $2";
			algebraic_solver=Reduce_Solver;			
			Solve='-A';
			shift;
			;;
		    *)
			echo $1 is an unknown solver - use hybrd, hooke or reduce;
			exit;;
	        esac;;
	-s )
                mtt_switches="$mtt_switches $1";
		sensitivity=sensitivity ;;
	-ss )
                mtt_switches="$mtt_switches $1";
		steadystate_computation=yes ;;
	-d )
1286
1287
1288
1289
1290
1291
1292

1293
1294
1295
1296
1297
1298
1299
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331







+







    echo '         -I  prints more information'
    echo '         -abg start at abg.m representation'
    echo '         -c  c-code generation'
    echo '         -d  <dir>  use directory <dir>'
    echo '         -dc Maximise derivative (not integral) causality'
    echo '         -dc Maximise derivative (not integral) causality'
    echo '         -i <implicit|euler|rk4>   Use implicit, euler or rk4 integration'
    echo '         -ae <reduce|hybrd|hooke>   Solve algebraic equations with Reduce, hybrd (fsolve) or "Hooke and Jeeves"'
    echo '         -o ode is same as dae'
    echo '         -oct use oct files in place of m files where appropriate'
    echo '         -opt optimise code generation'
    echo '         -p  print environment variables'
    echo '         -partition partition hierachical system'
    echo '         -pdf generate pdf in place of ps'
    echo '         -r  reset time stamp on representation'
1923
1924
1925
1926
1927
1928
1929



1930
1931
1932
1933
1934
1935
1936
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971







+
+
+







.PRECIOUS: %.cc # Don't let mtt delete them
$1_%.cc:  $1_%.m
	mtt_m2cc.sh  $1 \$* cat 

mtt_%.cc: ${MTT_LIB}/cc/mtt_%.cc
	cp ${MTT_LIB}/cc/mtt_\$*.cc mtt_\$*.cc

mtt_%.hh: ${MTT_LIB}/cc/mtt_%.hh
	cp ${MTT_LIB}/cc/mtt_\$*.hh mtt_\$*.hh

## .o files
.PRECIOUS: $1_%.o
$1_%.o: $1_%.cc $1_def.h $1_sympar.h $1_cr.h
	echo Compiling $1_\$*.cc
	${MTT_CXX} ${MTT_CXXFLAGS} ${MTT_CXXINCS} -c $< -DSTANDALONE

.PRECIOUS: mtt_%.o
2225
2226
2227
2228
2229
2230
2231





2232
2233
2234
2235
2236
2237
2238
2260
2261
2262
2263
2264
2265
2266
2267
2268
2269
2270
2271
2272
2273
2274
2275
2276
2277
2278







+
+
+
+
+







	@echo > /dev/null
$1_ode2odes_implicit.o  : $1_cseo.o    $1_csex.o   $1_smxa.o   $1_smxax.o   mtt_implicit.o
	@echo "Creating $1_ode2odes_implicit.o"
	ar -cr \$@ \$^
$1_ode2odes_implicit.oct: $1_cseo.oct  $1_csex.oct $1_smxa.oct $1_smxax.oct mtt_implicit.oct
	@echo > /dev/null

mtt_Solver.cc:	mtt_Solver.hh
$1_ode2odes_${algebraic_solver}.cc:	mtt_Solver.cc mtt_${algebraic_solver}.hh mtt_${algebraic_solver}.cc
$1_ode2odes_${algebraic_solver}.o:	mtt_Solver.o mtt_${algebraic_solver}.o
	@echo "Creating $1_ode2odes_${algebraic_solver}.o"
	ar -cr \$@ \$^

#SUMMARY numpar	numerical parameter declaration (m) 
$1_numpar.m:  $1_numpar.txt $1_sympars.txt
	mtt_txt2m $1 numpar

#SUMMARY numpar	numerical parameter declaration (c) 
#SUMMARY numpar	numerical parameter declaration (view) 
2527
2528
2529
2530
2531
2532
2533
2534

2535
2536
2537
2538
2539
2540
2541
2542
2543
2544
2545
2546

2547
2548
2549
2550
2551
2552
2553
2554
2555
2556

2557
2558
2559

2560
2561
2562

2563
2564
2565
2566
2567
2568
2569
2570
2571
2572
2573

2574
2575

2576
2577
2578
2579
2580
2581
2582
2567
2568
2569
2570
2571
2572
2573

2574
2575
2576
2577
2578
2579
2580
2581
2582
2583
2584
2585

2586
2587
2588
2589
2590
2591
2592
2593
2594
2595

2596
2597
2598

2599
2600
2601

2602
2603
2604
2605
2606
2607
2608
2609
2610
2611
2612

2613
2614

2615
2616
2617
2618
2619
2620
2621
2622







-
+











-
+









-
+


-
+


-
+










-
+

-
+







                $1_smxa.m $1_smxax.m\
                $1_simpar.m $1_numpar.m $1_state.m $1_input.m \
                $1_csex.m $1_cseo.m  $1_logic.m
ifeq ($using_oct,yes)
	touch $1_ode2odes.m # Create a dummy which wont' be used
	mtt $mtt_switches -q -u $1 ode2odes oct
else
	make_ode2odes $1 m $integration_method
	make_ode2odes $1 m $integration_method $algebraic_solver
endif
endif
ifneq ($integration_method,implicit)
$1_ode2odes.m : $1_def.r $1_sympars.txt\
		$1_simpar.m $1_numpar.m $1_state.m $1_input.m \
		$1_ode.m $1_odeo.m  $1_logic.m
ifeq ($using_oct,yes)
	echo "*** Warning: Shouldn't be here! Creating dummy $1_ode2odes.m"
	touch $1_ode2odes.m # Create a dummy which wont' be used
	mtt $mtt_switches -q -u $1 ode2odes oct
else
	make_ode2odes $1 m $integration_method
	make_ode2odes $1 m $integration_method $algebraic_solver
endif

endif

#SUMMARY ode2odes Simulation function (m)
#SUMMARY ode2odes Simulation function (cc)
#SUMMARY ode2odes Simulation function (oct)
#SUMMARY ode2odes Simulation function (exe)
$1_ode2odes.exe: $1_def.h $1_sympar.h\
		 $1_ode2odes.o $1_ode2odes_common.o $1_ode2odes_${integration_method}.o
		 $1_ode2odes.o $1_ode2odes_common.o $1_ode2odes_${integration_method}.o $1_ode2odes_${algebraic_solver}.o
	echo Creating $1_ode2odes.exe
	${MTT_CXX} ${MTT_CXXFLAGS} -o $1_ode2odes.exe\
	$1_ode2odes.o $1_ode2odes_common.o $1_ode2odes_${integration_method}.o\
	$1_ode2odes.o $1_ode2odes_common.o $1_ode2odes_${integration_method}.o $1_ode2odes_${algebraic_solver}.o\
	${MTT_LDFLAGS} ${MTT_CXXLIBS}

$1_ode2odes.o: $1_ode2odes.cc $1_ode2odes_common.o $1_ode2odes_${integration_method}.o
$1_ode2odes.o: $1_ode2odes.cc $1_ode2odes_common.o $1_ode2odes_${integration_method}.o $1_ode2odes_${algebraic_solver}.o
	echo Creating $1_ode2odes.o
	${MTT_CXX} ${MTT_CXXFLAGS} ${MTT_CXXINCS} -c $1_ode2odes.cc -DSTANDALONE

$1_ode2odes.oct: $1_ode2odes.cc $1_ode2odes_common.oct $1_ode2odes_${integration_method}.oct
	touch $1_ode2odes.m
	echo Creating $1_ode2odes.oct
	$MKOCTFILE $1_ode2odes.cc

$1_ode2odes.cc: $1_def.r $1_sympars.txt\
		$1_ode2odes_common.m $1_ode2odes_common.cc\
		$1_ode2odes_${integration_method}.m $1_ode2odes_${integration_method}.cc
		$1_ode2odes_${integration_method}.m $1_ode2odes_${integration_method}.cc mtt_Solver.cc mtt_${algebraic_solver}.cc mtt_${algebraic_solver}.hh
	touch $1_ode2odes.m
	make_ode2odes $1 cc $integration_method
	make_ode2odes $1 cc $integration_method $algebraic_solver

#Conversion of m to p to c
#SUMMARY ode	ordinary differential equations (c)
#SUMMARY ode	ordinary differential equations (p)
#SUMMARY state	state declaration (c) 
#SUMMARY state	state declaration (p) 
$1_simpar.p : $1_def.r $1_simpar.m

Modified mttroot/mtt/bin/trans/make_ode2odes from [afde7a9c2a] to [a0a317a4de].

1
2
3
4
5
6
7
8
9
10
11





12
13
14
15
16
17
18
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23











+
+
+
+
+







#! /bin/sh

     ###################################### 
     ##### Model Transformation Tools #####
    ######################################

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.57.2.1  2001/05/04 04:07:24  geraint
## Numerical solution of algebraic equations.
## sys_ae.cc written for unsolved inputs.
## Solution of equations using hybrd from MINPACK (as used by Octave fsolve).
##
## Revision 1.57  2001/04/01 03:38:54  geraint
## Reset row to zero after write to file, ready for subsequent runs.
## Eliminates SIGSEGV in Octave when _ode2odes called multiple times.
##
## Revision 1.56  2001/03/30 15:13:58  gawthrop
## Rationalised simulation modes to each return mtt_data
##
221
222
223
224
225
226
227






228
229
230
231
232
233
234
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245







+
+
+
+
+
+







filename=${sys}_ode2odes.${lang}

if [ -n "$3" ]; then
  method=$3    
else
  method=implicit  
fi

if [ -n "$4" ]; then
    algebraic_solver=$4
else
  algebraic_solver="Reduce_Solver"
fi

echo Creating $filename with $method integration method

# Find system constants
Nx=`mtt_getsize $sys x` # States
Nu=`mtt_getsize $sys u` # Inputs 
Ny=`mtt_getsize $sys y` # Inputs  
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366


367
368
369
370
371
372
373
357
358
359
360
361
362
363

364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385







-













+
+







	;;
esac

cat <<EOF  > $filename
#include <octave/oct.h>
#include <octave/load-save.h>
#include <octave/lo-mappers.h>
#include <octave/NLEqn.h>
#include <octave/ov-struct.h>
#include <octave/variables.h>

#ifndef STANDALONE
#include <octave/${feval_header}>
#endif

#include "${sys}_def.h"
#include "${sys}_sympar.h"

#ifdef STANDALONE
#include <csignal>
#include <fstream>

#include "mtt_${algebraic_solver}.hh"

extern ColumnVector F${sys}_ae (
	ColumnVector &x,
	ColumnVector &u,
	const double &t,
	ColumnVector &par);

442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557


558
559
560
561
562
563
564
454
455
456
457
458
459
460


































































































461
462
463
464
465
466
467
468
469
470

471
472
473
474
475
476
477
478
479







-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-










-
+
+







fi
cat <<EOF >> $filename

void set_signal_handlers (void);

#endif // STANDALONE

class Solver
{
public:
  Solver (void) { 
    solve_errors.input		= 0;
    solve_errors.user		= 0;
    solve_errors.converge	= 0;
    solve_errors.progress	= 0;
    solve_errors.limit		= 0;
    solve_errors.unknown	= 0;
 }
  ColumnVector solve (const ColumnVector &x,
		      const ColumnVector &u,
		      const double &t,
		      const ColumnVector &par)
  {
    initialise(x,u,t,par);
    NLFunc fcn(&Solver::evaluate);
    NLEqn  eqn(Ui,fcn);
    Ui = eqn.solve(info);
    switch (info)
    {
      case -2:
	solve_errors.input++;
	write_stats();
        break;
      case -1:
	solve_errors.user++;
	write_stats();
        break;
      case 1:
	solve_errors.converge++;
        break;
      case 3:
	solve_errors.progress++;
	write_stats();
        break;
      case 4:
	solve_errors.limit++;
	write_stats();
	break;
      default:
	solve_errors.unknown++;
	write_stats();
	break;;
    }
    U.insert (Ui,MTTNU);
    return U;
  }
  void write_stats (void)
  {
    cerr << "input error (" << solve_errors.input << ") "
	 << ",  user (" << solve_errors.user << ") "
	 << ",  converge (" << solve_errors.converge << ") "
	 << ",  progress (" << solve_errors.progress << ") "
	 << ",  limit (" << solve_errors.limit << ") "
	 << ",  unknown (" << solve_errors.unknown << ") "
	 << endl;
  }
private:
  static ColumnVector evaluate (const ColumnVector &tryUi)
  {
    U.insert(tryUi,MTTNU);
    return F${sys}_ae(X,U,T,P);
  }
  void initialise (const ColumnVector &x,
		   const ColumnVector &u,
		   const double &t,
		   const ColumnVector &par)
  {
    X = x;
    T = t;
    P = par;
    U.insert(u,0);
    Ui = ColumnVector(MTTNYZ,0.0);
  }
private:
  struct {
    long int input;
    long int user;
    long int converge;
    long int progress;
    long int limit;
    long int unknown; 
  } solve_errors;
  static double T;
  static ColumnVector X;
  static ColumnVector P;
  static ColumnVector U;
  static ColumnVector Ui;
  int info;
};
double Solver::T = 0.0;
ColumnVector Solver::U(MTTNU+MTTNYZ,0.0);
ColumnVector Solver::X(MTTNX,0.0);
ColumnVector Solver::P(MTTNPAR,0.0);
ColumnVector Solver::Ui(MTTNYZ,0.0);

inline ColumnVector
mtt_input (ColumnVector &x,
	   ColumnVector &y,
	   const double &t,
	   ColumnVector &par)
{
#ifdef STANDALONE
  static ColumnVector u  (MTTNU);
  static ColumnVector ui (MTTNYZ);
  static ColumnVector U  (MTTNU+MTTNYZ);
  static Solver ae;
 
  static ${algebraic_solver} ae(F${sys}_ae,MTTNPAR,MTTNU,MTTNX,MTTNY,MTTNYZ);

  u = F${sys}_input (x, y, t, par);
  if (MTTNYZ == 0)
    {
      return u;
    }
  else

Modified mttroot/mtt/cc/include/useful-functions.hh from [61119a71b8] to [652bf845ed].

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








-
-
-




-
+





-
+





-
+





-
+







#define pi 3.14159 // Predefine pi

#ifndef HAVE_USEFUL_FUNCTIONS_HH
#define HAVE_USEFUL_FUNCTIONS_HH


#ifndef __cplusplus
#define inline			/* strip */
typedef double  doubleref_t;
#else
typedef double &doubleref_t;
#endif // ! __cplusplus


static inline double
max (const doubleref_t x1, const doubleref_t x2)
max (const double x1, const double x2)
{
  return static_cast<double>((x1 >= x2) ? x1 : (x1 < x2) ? x2 : 0);
}

static inline double
min (const doubleref_t x1, const doubleref_t x2)
min (const double x1, const double x2)
{
  return static_cast<double>((x1 <= x2) ? x1 : (x1 > x2) ? x2 : 0);
}

static inline double
nonsingular (const doubleref_t x)
nonsingular (const double x)
{
  return static_cast<double>((x == 0) ? 1.0e-30 : x);
}

static inline double
sign (const doubleref_t x)
sign (const double x)
{
  return static_cast<double>((x > 0) ? +1 : (x < 0) ? -1 : 0);
}


// Octave functions
#ifdef __cplusplus

Added mttroot/mtt/lib/cc/mtt_HJ_Solver.cc version [649cd58504].
























































































































































































































































































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
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
113
114
115
116
117
118
119
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
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
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

#include "mtt_HJ_Solver.hh"

// http://www.netlib.org/opt/hooke.c
// Hooke and Jeeves solution
  
HJ_Solver *HJ_Solver::static_ptr;

double
HJ_Solver::f (double tryUi[], int nyz)
{
  for (int i = 0; i < nyz; i++)
    HJ_Solver::static_ptr->_ui(i) = tryUi[i];

  HJ_Solver::static_ptr->_yz = HJ_Solver::static_ptr->eval(HJ_Solver::static_ptr->_ui);

  double ans = 0.0;
  for (int i = 0; i < nyz; i++)
    ans += HJ_Solver::static_ptr->_yz(i)*HJ_Solver::static_ptr->_yz(i);
  return ans;
}

void
HJ_Solver::Solve (void)
{
  double startpt[_nyz];
  double endpt[_nyz];
  for (int i = 0; i < _nyz; i++)
    startpt[i] = _ui(i);
  hooke(_nyz,startpt,endpt);
  for (int i = 0; i < _nyz; i++)
    _ui(i) = endpt[i];
}

  // adapted from http://www.netlib.org/opt/hooke.c:

/* Nonlinear Optimization using the algorithm of Hooke and Jeeves  */
/*      12 February 1994        author: Mark G. Johnson            */


/* Find a point X where the nonlinear function f(X) has a local    */
/* minimum.  X is an n-vector and f(X) is a scalar.  In mathe-     */
/* matical notation  f: R^n -> R^1.  The objective function f()    */
/* is not required to be continuous.  Nor does f() need to be      */
/* differentiable.  The program does not use or require            */
/* derivatives of f().                                             */

/* The software user supplies three things: a subroutine that      */
/* computes f(X), an initial "starting guess" of the minimum point */
/* X, and values for the algorithm convergence parameters.  Then   */
/* the program searches for a local minimum, beginning from the    */
/* starting guess, using the Direct Search algorithm of Hooke and  */
/* Jeeves.                                                         */

/* This C program is adapted from the Algol pseudocode found in    */
/* "Algorithm 178: Direct Search" by Arthur F. Kaupe Jr., Commun-  */
/* ications of the ACM, Vol 6. p.313 (June 1963).  It includes the */
/* improvements suggested by Bell and Pike (CACM v.9, p. 684, Sept */
/* 1966) and those of Tomlin and Smith, "Remark on Algorithm 178"  */
/* (CACM v.12).  The original paper, which I don't recommend as    */
/* highly as the one by A. Kaupe, is:  R. Hooke and T. A. Jeeves,  */
/* "Direct Search Solution of Numerical and Statistical Problems", */
/* Journal of the ACM, Vol. 8, April 1961, pp. 212-229.            */

/* Calling sequence:                                               */
/*  int hooke(nvars, startpt, endpt, rho, epsilon, itermax)        */
/*                                                                 */
/*     nvars       {an integer}  This is the number of dimensions  */
/*                 in the domain of f().  It is the number of      */
/*                 coordinates of the starting point (and the      */
/*                 minimum point.)                                 */
/*     startpt     {an array of doubles}  This is the user-        */
/*                 supplied guess at the minimum.                  */
/*     endpt       {an array of doubles}  This is the location of  */
/*                 the local minimum, calculated by the program    */
/*     rho         {a double}  This is a user-supplied convergence */
/*                 parameter (more detail below), which should be  */
/*                 set to a value between 0.0 and 1.0.  Larger     */
/*                 values of rho give greater probability of       */
/*                 convergence on highly nonlinear functions, at a */
/*                 cost of more function evaluations.  Smaller     */
/*                 values of rho reduces the number of evaluations */
/*                 (and the program running time), but increases   */
/*                 the risk of nonconvergence.  See below.         */
/*     epsilon     {a double}  This is the criterion for halting   */
/*                 the search for a minimum.  When the algorithm   */
/*                 begins to make less and less progress on each   */
/*                 iteration, it checks the halting criterion: if  */
/*                 the stepsize is below epsilon, terminate the    */
/*                 iteration and return the current best estimate  */
/*                 of the minimum.  Larger values of epsilon (such */
/*                 as 1.0e-4) give quicker running time, but a     */
/*                 less accurate estimate of the minimum.  Smaller */
/*                 values of epsilon (such as 1.0e-7) give longer  */
/*                 running time, but a more accurate estimate of   */
/*                 the minimum.                                    */
/*     itermax     {an integer}  A second, rarely used, halting    */
/*                 criterion.  If the algorithm uses >= itermax    */
/*                 iterations, halt.                               */


/* The user-supplied objective function f(x,n) should return a C   */
/* "double".  Its  arguments are  x -- an array of doubles, and    */
/* n -- an integer.  x is the point at which f(x) should be        */
/* evaluated, and n is the number of coordinates of x.  That is,   */
/* n is the number of coefficients being fitted.                   */

/* rho, the algorithm convergence control                          */
/*      The algorithm works by taking "steps" from one estimate of */
/*    a minimum, to another (hopefully better) estimate.  Taking   */
/*    big steps gets to the minimum more quickly, at the risk of   */
/*    "stepping right over" an excellent point.  The stepsize is   */
/*    controlled by a user supplied parameter called rho.  At each */
/*    iteration, the stepsize is multiplied by rho  (0 < rho < 1), */
/*    so the stepsize is successively reduced.                     */
/*      Small values of rho correspond to big stepsize changes,    */
/*    which make the algorithm run more quickly.  However, there   */
/*    is a chance (especially with highly nonlinear functions)     */
/*    that these big changes will accidentally overlook a          */
/*    promising search vector, leading to nonconvergence.          */
/*      Large values of rho correspond to small stepsize changes,  */
/*    which force the algorithm to carefully examine nearby points */
/*    instead of optimistically forging ahead.  This improves the  */
/*    probability of convergence.                                  */
/*      The stepsize is reduced until it is equal to (or smaller   */
/*    than) epsilon.  So the number of iterations performed by     */
/*    Hooke-Jeeves is determined by rho and epsilon:               */
/*          rho**(number_of_iterations) = epsilon                  */
/*      In general it is a good idea to set rho to an aggressively */
/*    small value like 0.5 (hoping for fast convergence).  Then,   */

  
/*    if the user suspects that the reported minimum is incorrect  */
/*    (or perhaps not accurate enough), the program can be run     */
/*    again with a larger value of rho such as 0.85, using the     */
/*    result of the first minimization as the starting guess to    */
/*    begin the second minimization.                               */

/* Normal use: (1) Code your function f() in the C language        */
/*             (2) Install your starting guess {or read it in}     */
/*             (3) Run the program                                 */
/*             (4) {for the skeptical}: Use the computed minimum   */
/*                    as the starting point for another run        */

/* Data Fitting:                                                   */
/*      Code your function f() to be the sum of the squares of the */
/*      errors (differences) between the computed values and the   */
/*      measured values.  Then minimize f() using Hooke-Jeeves.    */
/*      EXAMPLE: you have 20 datapoints (ti, yi) and you want to   */
/*      find A,B,C such that  (A*t*t) + (B*exp(t)) + (C*tan(t))    */
/*      fits the data as closely as possible.  Then f() is just    */
/*      f(x) = SUM (measured_y[i] - ((A*t[i]*t[i]) + (B*exp(t[i])) */
/*                                + (C*tan(t[i]))))^2              */
/*      where x[] is a 3-vector consisting of {A, B, C}.           */

/*                                                                 */
/*  The author of this software is M.G. Johnson.                   */
/*  Permission to use, copy, modify, and distribute this software  */
/*  for any purpose without fee is hereby granted, provided that   */
/*  this entire notice is included in all copies of any software   */
/*  which is or includes a copy or modification of this software   */
/*  and in all copies of the supporting documentation for such     */
/*  software.  THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT    */
/*  ANY EXPRESS OR IMPLIED WARRANTY.  IN PARTICULAR, NEITHER THE   */
/*  AUTHOR NOR AT&T MAKE ANY REPRESENTATION OR WARRANTY OF ANY     */
/*  KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS    */
/*  FITNESS FOR ANY PARTICULAR PURPOSE.                            */
/*                                                                 */

double
HJ_Solver::best_nearby(double *delta,
		       double *point,
		       double prevbest,
		       int    nvars)
{
  double          z[VARS];
  double          minf, ftmp;
  int             i;
  minf = prevbest;
  for (i = 0; i < nvars; i++)
    z[i] = point[i];
  for (i = 0; i < nvars; i++) {
    z[i] = point[i] + delta[i];
    ftmp = f(z, nvars);
    if (ftmp < minf)
      minf = ftmp;
    else {
      delta[i] = 0.0 - delta[i];
      z[i] = point[i] + delta[i];
      ftmp = f(z, nvars);
      if (ftmp < minf)
	minf = ftmp;
      else
	z[i] = point[i];
    }
  }
  for (i = 0; i < nvars; i++)
    point[i] = z[i];
  return (minf);
}


int
HJ_Solver::hooke (int	nvars,			// MTTNYZ
		  double	startpt[],	// user's initial guess
		  double	endpt[],	// result
		  double	rho,		// geometric shrink factor
		  double	epsilon,	// end value stepsize
		  int	itermax)		// max # iterations 
{
  const int VARS = _nyz;
  double          delta[VARS];
  double          newf, fbefore, steplength, tmp;
  double          xbefore[VARS], newx[VARS];
  int             i, j, keep;
  int             iters, iadj;
  for (i = 0; i < nvars; i++) {
    newx[i] = xbefore[i] = startpt[i];
    delta[i] = fabs(startpt[i] * rho);
    if (delta[i] == 0.0)
      delta[i] = rho;
  }
  iadj = 0;
  steplength = rho;
  iters = 0;
  fbefore = f(newx, nvars);
  newf = fbefore;
  while ((iters < itermax) && (steplength > epsilon)) {
    iters++;
    iadj++;
  /* find best new point, one coord at a time */
  for (i = 0; i < nvars; i++) {
    newx[i] = xbefore[i];
  }
  newf = best_nearby(delta, newx, fbefore, nvars);
  /* if we made some improvements, pursue that direction */
  keep = 1;
  while ((newf < fbefore) && (keep == 1)) {
    iadj = 0;
    for (i = 0; i < nvars; i++) {
      /* firstly, arrange the sign of delta[] */
      if (newx[i] <= xbefore[i])
	delta[i] = 0.0 - fabs(delta[i]);
      else
	delta[i] = fabs(delta[i]);                                   /* now, move further in this direction */
      tmp = xbefore[i];
      xbefore[i] = newx[i];
      newx[i] = newx[i] + newx[i] - tmp;
    }
    fbefore = newf;
    newf = best_nearby(delta, newx, fbefore, nvars);
    /* if the further (optimistic) move was bad.... */
    if (newf >= fbefore)
      break;
    /* make sure that the differences between the new */
    /* and the old points are due to actual */
    /* displacements; beware of roundoff errors that */
    /* might cause newf < fbefore */
    keep = 0;
    for (i = 0; i < nvars; i++) {
      keep = 1;
      if (fabs(newx[i] - xbefore[i]) >
	  (0.5 * fabs(delta[i])))
	break;
      else
	keep = 0;
    }
  }
  if ((steplength >= epsilon) && (newf >= fbefore)) {
    steplength = steplength * rho;
    for (i = 0; i < nvars; i++) {
      delta[i] *= rho;
    }
  }
}
for (i = 0; i < nvars; i++)
  endpt[i] = xbefore[i];
return (iters);
}

Added mttroot/mtt/lib/cc/mtt_HJ_Solver.hh version [da2fe34a04].




















































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

#include "mtt_Solver.hh"

class HJ_Solver : public Solver {

  // http://www.netlib.org/opt/hooke.c
  // Hooke and Jeeves solution
  
public:
  
  HJ_Solver (sys_ae ae,
	     const int npar,
	     const int nu,
	     const int nx,
	     const int ny,
	     const int nyz)
    : Solver (ae,npar,nu,nx,ny,nyz)
  { static_ptr = this; VARS = nyz; };
  
  static double
  f (double tryUi[], int nyz);
  
protected:

  void
  Solve (void);

  double
  best_nearby (double    	delta[],
	       double   	point[],
	       double   	prevbest,
	       int      	nvars);
  
  int
  hooke (int	nvars,			  // MTTNYZ
	 double	startpt[],		  // user's initial guess
	 double	endpt[],		  // result
	 double	rho		= 0.05,	  // geometric shrink factor
	 double	epsilon		= 1e-3,	  // end value stepsize
	 int	itermax		= 5000);  // max # iterations 

private:

  int VARS;

public:

  static HJ_Solver *static_ptr;

};

Added mttroot/mtt/lib/cc/mtt_Hybrd_Solver.cc version [e965430e5f].




































































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
60
61
62
63
64
65
66
67
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

#include "mtt_Hybrd_Solver.hh"

// http://www.netlib.org/minpack/hybrd.f
// used by Octave's fsolve
  
Hybrd_Solver *Hybrd_Solver::static_ptr;

ColumnVector
Hybrd_Solver::f_hybrd (const ColumnVector &tryUi)
{
  Hybrd_Solver::static_ptr->_yz = Hybrd_Solver::static_ptr->eval(tryUi);
  return Hybrd_Solver::static_ptr->_yz;
}

void
Hybrd_Solver::Solve (void)
{    
  int info;
  static int input_errors;
  static int user_errors;
  static int convergences;
  static int progress_errors;
  static int limit_errors;
  static int unknown_errors;
  
  NLFunc fcn(&Hybrd_Solver::f_hybrd);
  NLEqn	 eqn(Solver::_ui,fcn);
  //  eqn.set_tolerance(0.01);
  Solver::_ui = eqn.solve(info);

  switch (info)
    {
    case 1:
      convergences++;
      break;
    case -2:
      input_errors++;
      break;
    case -1:
      user_errors++;
      break;
    case 3:
      progress_errors++;
      break;
    case 4:
      //      if (abs(eval(_ui).max()) > 1.0e-6)
      limit_errors++;
      //      else
      //	convergences++;
      break;
    default:
      unknown_errors++;
      break;
    }
  if (1 != info)
    {
      cerr << "input (" << input_errors << ") "
	   << "  user (" << user_errors << ") "
	   << "  converge (" << convergences << ") "
	   << "  progress (" << progress_errors << ") "
	   << "  limit (" << limit_errors << ")"
	   << "  unknown (" << unknown_errors << ")"
	   << "  (max error = " << abs(eval(_ui).max()) << ")" << endl;
    }
}

Added mttroot/mtt/lib/cc/mtt_Hybrd_Solver.hh version [79adf97ac2].





































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

#include "mtt_Solver.hh"
#include <octave/NLEqn.h>

class Hybrd_Solver : public Solver {

  // http://www.netlib.org/minpack/hybrd.f
  // used by Octave's fsolve
  
public:

  Hybrd_Solver (sys_ae ae,
		const int npar,
		const int nu,
		const int nx,
		const int ny,
		const int nyz)
  : Solver (ae,npar,nu,nx,ny,nyz)
  {
    static_ptr = this;
  }

  static ColumnVector
  f_hybrd (const ColumnVector &tryUi);

protected:

  void
  Solve (void);

public:

  static Hybrd_Solver *static_ptr;

};

Added mttroot/mtt/lib/cc/mtt_Reduce_Solver.cc version [55dcc109c3].












1
2
3
4
5
6
7
8
9
10
11
+
+
+
+
+
+
+
+
+
+
+

#include "mtt_Reduce_Solver.hh"


void
Reduce_Solver::Solve (void)
{
  cerr << "Error:"
       << " Symbolic solution of equations failed during model build" << endl
       << "       Try using one of the other algebraic solution methods" << endl;
}

Added mttroot/mtt/lib/cc/mtt_Reduce_Solver.hh version [b44a93cb60].


























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

#include "mtt_Solver.hh"

class Reduce_Solver : public Solver {

  // Dummy class
  // This will not be used unless the Reduce solver has failed earlier
  // in the model build process

public:

  Reduce_Solver (sys_ae ae,
		 const int npar,
		 const int nu,
		 const int nx,
		 const int ny,
		 const int nyz)
    : Solver (ae,npar,nu,nx,ny,nyz)
  { ; };
	
  void
  Solve (void);

};
   

Added mttroot/mtt/lib/cc/mtt_Solver.cc version [d33ce07dea].











































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

#include "mtt_Solver.hh"

Solver::Solver (sys_ae ae,
		const int npar,
		const int nu,
		const int nx,
		const int ny,
		const int nyz)
{
  _ae  = ae;
  _np  = npar; 
  _nu  = nu;
  _nx  = nx;
  _ny  = ny;
  _nyz = nyz;

  _uui = ColumnVector (_nu+_nyz);
};

ColumnVector
Solver::solve (const ColumnVector	&x,
	       const ColumnVector	&u,
	       const double		&t,
	       const ColumnVector	&par)
{
  _x = x;
  _uui.insert(u,0);
  _t = t;
  _par = par;
  _ui = ColumnVector(_nyz,0.2);
  Solve ();
  _uui.insert(_ui,_nu);
  return _uui;
}    

ColumnVector
Solver::eval (const ColumnVector &ui)
{
  _uui.insert(ui,_nu);
  return _ae (_x, _uui, _t, _par);
}

Added mttroot/mtt/lib/cc/mtt_Solver.hh version [d701ed61c9].





























































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

#ifndef MTT_SOLVER
#define MTT_SOLVER

#include <cmath>
#include <cstdlib>
#include <iostream>

#include <octave/oct.h>

class Solver {

  typedef ColumnVector (*sys_ae) // pointer to F${sys}_ae function
    (ColumnVector &,ColumnVector &,const double &t,ColumnVector &);

public:

  Solver (sys_ae ae,
	  const int npar,
	  const int nu,
	  const int nx,
	  const int ny,
	  const int nyz);

  ColumnVector
  solve (const ColumnVector	&x,
	 const ColumnVector	&u,
	 const double		&t,
	 const ColumnVector	&par);
  
  ColumnVector
  eval (const ColumnVector	&ui);
  
protected:
  
  virtual void
  Solve (void) = 0;

protected:

  ColumnVector			_x;
  ColumnVector			_uui;
  double			_t;
  ColumnVector			_par;

  ColumnVector  		_ui;
  ColumnVector          	_yz;

  int   			_nu;
  int				_np;
  int				_nx;
  int				_ny;
  int				_nyz;
  
  sys_ae			_ae;

};

#endif // MTT_SOLVER


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