Optimization: principles and algorithms, by Michel Bierlaire
torczon.m
Go to the documentation of this file.
1 %> \file
2 %> Algorithm 15.2: Torczon algorithm. Implementation of algorithm 15.2 of \cite Bier15-book
3 %>
4 %> @note Tested with \ref run1503torczon.m
5 %>
6 %> @author Michel Bierlaire
7 %> @date Sun Mar 22 10:26:50 2015
8 %> @ingroup Algorithms
9 %> @ingroup chap15
10 
11 %> Applies Torczon's algorithm to solve \f$\min_x f(x)\f$ where \f$f:\mathbb{R}^n\to\mathbb{R}\f$.
12 %> @param obj the name of the Octave function defining f(x)
13 %> @param x0 matrix nx(n+1) containnig the initial simplex
14 %> @param eps algorithm stops if \f$\|x_{n+1}-x_1\| \leq \varepsilon \f$.
15 %> @param maxiter maximum number of iterations (Default: 100)
16 %> @return [solution,simplex]: the solution and the final simplex.
17 function [solution,simplex] = torczon(obj,x0,eps,maxiter=100)
18  x = x0 ;
19  [n,m] = size(x0) ;
20  if (m != n+1)
21  error("Incorrect size") ;
22  endif
23  iter = 0 ;
24  do
25  iter = iter + 1 ;
26  f = [] ;
27  fr = [] ;
28  fe = [] ;
29  for k = 1:n+1
30  f = [ f ; feval(obj,x(:,k)) ] ;
31  endfor
32 
33  ## Worst value
34  [worst,worstindex] = max(f) ;
35  ## Best value
36  [best,bestindex] = min(f) ;
37 
38  for jj=1:n
39  printf("%e\t",x(jj,bestindex)) ;
40  endfor
41  printf("%e\t",best) ;
42 
43  d = x(:,worstindex) - x(:,bestindex) ;
44 
45  ## Reflexion
46  for k = 1:n+1
47  if (k == bestindex)
48  xr(:,k) = x(:,k) ;
49  else
50  xr(:,k) = 2 * x(:,bestindex) - x(:,k) ;
51  endif
52  endfor
53 
54  for k = 1:n+1
55  fr = [ fr ; feval(obj,xr(:,k)) ] ;
56  endfor
57 
58  ## Best value
59  [bestr,bestrindex] = min(fr) ;
60 
61  if (bestr >= best)
62  ## Contraction
63  printf("C\n") ;
64  for k = 1:n+1
65  if (k != bestindex)
66  x(:,k) = 0.5 * (x(:,bestindex) + x(:,k)) ;
67  endif
68  endfor
69  else
70  ## Expansion
71  for k = 1:n+1
72  if (k == bestindex)
73  xe(:,k) = x(:,k) ;
74  else
75  xe(:,k) = 3 * x(:,bestindex) - 2 * x(:,k) ;
76  endif
77  endfor
78 
79  for k = 1:n+1
80  fe = [ fe ; feval(obj,xe(:,k)) ] ;
81  endfor
82 
83  ## Best value
84  [beste,besteindex] = min(fe) ;
85 
86  if (beste >= bestr)
87  printf("R\n") ;
88  x = xr ;
89  else
90  printf("E\n") ;
91  x = xe ;
92  endif
93  endif
94 
95  until (norm(d) < eps || iter >= maxiter)
96  solution = x(:,bestindex) ;
97  simplex = x ;
98 endfunction
Copyright 2015-2016 Michel Bierlaire