Bài giảng chương 8: Tối ưu hóa

370 CHƯƠNG 8: TỐI ƯU HOÁ §1. KHÁI NIỆM CHUNG VỀ TỐI ƯU HOÁ    Tối ưu hoá là thuật ngữ thường được dùng để cực tiểu hoá hay cực đại  hoá một hàm. Thông thường ta chỉ cần tìm cực tiểu một hàm là đủ. Việc tìm  cực đại của f(x) thực hiện một cách đơn giản bằng cách tìm cực tiểu của hàm  −f(x) . Hàm f là hàm giá trị hay hàm đối tượng, cần được giữ cực tiểu. Biến x  là biến có thể hiệu chỉnh tự do.    Các thuật toán cực tiểu hoá là các thủ thuật lặp đòi hỏi một giá trị ban  đầu của biến x. Nếu f(x) có nhiều cực tiểu địa phương, việc chọn giá trị đầu sẽ  xác định cực tiểu nào được tính. Ta không có cách nào bảo đảm là tìm được  cực tiểu toàn cục.    Các biến  có  thể bị  ràng buộc bằng  các  đẳng  thức hay bất  đẳng  thức.  Phần lớn các phương pháp là tìm cực tiểu không ràng buộc, nghĩa là không có  hạn chế nào đối với biến x. Các bài toán này bao gồm tìm cực tiểu của hàm,  tìm điểm tĩnh ‐ điểm có gradient triệt tiêu. Các bài toán tìm cực tiểu có ràng  buộc khó hơn và thuật toán khá phức tạp.     Trong chương này chúng ta sẽ  lần  lượt xét các thuật toán tìm cực tiểu  không ràng buộc và có ràng buộc.  §2. PHƯƠNG PHÁP TIẾT DIỆN VÀNG   Ta xét bài toán tìm cực tiểu của hàm một biến f(x). Điểm cực tiểu được  xác định  theo điều kiện df/dx = 0. Do có  thể có nhiều điểm cực  tiểu nên  ta  phải vây điểm cực tiểu(xác định lân cận chứa điểm cực tiểu) trước. Thủ thuật  vây điểm cực tiểu khá đơn giản: cho điểm đầu x0 và tính giá trị của hàm đang  đi  xuống  tại  các  điểm  tiếp  theo  x1,  x2,  x3,...  cho  đến  tại  xn hàm  tăng  lại  thì  dừng. Điểm cực tiểu bị vây trong khoảng (xn‐2, xn). Khoảng (xi+1, xi) không nên  chọn là hằng số vì như vậy cần nhiều bước tính. Hợp lí nhất là nên tăng kích  thước bước tính để đạt được cực tiểu nhanh hơn, ngay cả khi cực tiểu bị vây  trong  một  đoạn  khá  rộng.  Ta  chọn  kích  thước  tăng  theo  dạng  hằng  số:  + =i 1 ih ch với  >c 1.  ta xây dựng hàm goldbracket() để vây điểm cực  tiểu của  hàm:  function [a, b] = goldbracket(func, x1, h)  % vay diem cuc tieu cua f(x).  c = 1.618033989;  371 f1 = feval(func, x1);  x2 = x1 + h;   f2 = feval(func,x2);  if f2 > f1      h = ‐h;      x2 = x1 + h;       f2 = feval(func, x2);      if f2 > f1          a = x2;           b = x1 ‐ h;           return      end  end  for i = 1:100      h = c*h;      x3 = x2 + h;       f3 = feval(func, x3);      if f3 > f2          a = x1;           b = x3;           return      end      x1 = x2;       f1 = f2;       x2 = x3;      f2 = f3;  end  error(ʹGoldbracket khong tim thay diem cuc tieuʹ)    Tiết diện vàng là một biến thể của phương pháp chia đôi dùng khi tìm  nghiệm của phương trình f(x) = 0. Giả sử điểm cực tiểu bị vây trong khoảng  (a,  b)  có  độ  dài  h.  Để  thu  nhỏ  khoảng  (a,  b)  ta  tính  giá  trị  của  hàm  tại  = −1x b rh  và  = +2x a rh  như hình vẽ. Nếu f1 = f(x1) lớn hơn f2 = f(x2) như hình  a  thì  cực  tiểu nằm  trong  khoảng  (x1,  b) nếu ngược  lại  cực  tiểu nằm  trong  khoảng (a, x2).  372   Giả sử f1 > f2 , ta đặt a = x1 và và x1 =  x2 và có khoảng (a, b) mới như hình b. Để  thực  hiện  bước  thu  gọn  tiếp  theo  ta  lại  tính giá trị của hàm tại x2 = a + rh’ và lặp  lại quá  trình. Quá  trình  làm việc chỉ nếu  hình a và hình b tương tự, nghĩa  là hằng  số r không đổi khi xác định x1 và x2 ở cả  hai hình. Từ hình a ta thấy:    − = −2 1x x 2rh h   Cùng một khoảng cách đó từ hình b ta có:    x1 ‐ a = h’ ‐ rh’  Cân bằng các khoảng này ta được:    2rh ‐ h = h’ ‐ rh’  Thay h’ = rh và khử h:    2r ‐ 1 = r(1 ‐ r)  Giải phương trình này ta nhận được tỉ lệ vàng:    −= =5 1r 0.618033989... 2 Chú ý là mỗi lần thu gọn khoảng chứa điểm cực tiểu thì khoảng (a, b) giảm tỉ  lệ với r. Điều này làm số lần tính lớn hơn phương pháp chia đôi. Tuy nhiên  phương pháp tỉ lệ vàng chỉ đòi hỏi tính giá trị hàm một lần trong khi phương  pháp chia đôi cần tính giá trị hàm 2 lần. Số lần tính xác định bằng:    − = εnb a r   hay:  ε − ε= = − − ln b an 2.078087n ln b a     h = b ‐ a = 0.382  Ta xây dựng hàm golden() để thực hiện thuật toán này:  function [xmin, ymin] = golden(f, a, b, delta, epsilon)  % a va b la doan tim cuc tieu  % delta sai so cua x  % epsilon sai so cua y  r1 = (sqrt(5) ‐ 1)/2;  r2 = r1^2;  h = b ‐ a;  a bx1  x2  h  rh  2rh‐h  rh  a  a bx1  x2  h’  rh’  rh’  b  373 fa = f(a);  fb = f(b);  x1 = a + r2*h;  x2 = a + r1*h;  f1 = f(x1);  f2 = f(x2);  k = 1;  while (abs(fb‐fa) > epsilon)|(h > delta)     k = k + 1;     if (f1 < f2)        b = x2;        fb = f2;        x2 = x1;        f2 = f1;        h = b ‐ a;        x1 = a + r2*h;        f1 = f(x1);     else        a = x1;        fa = f1;        x1 = x2;        f1 = f2;        h = b ‐ a;        x2 = a + r1*h;        f2 = f(x2);     end  end  dp = abs(b ‐ a);  dy = abs(fb ‐ fa);  p = a;  yp = fa;  if (fb < fa)     p = b;     yp = fb;  end  xmin = p;  374 ymin = yp;  Để tìm cực tiểu của hàm ta dùng chương trình ctgolden.m:  clear all, clc  f = inline(ʹ1.6*x^3 + 3*x^2 ‐2*xʹ);  x = 0;   delta = 1e‐8;  epsilon = 1e‐10;  [a, b] = goldbracket(f, x, 0.2);  [xmin, ymin] = golden(f, a, b, delta, epsilon)  §3. PHƯƠNG PHÁP XẤP XỈ BẬC HAI    Ý tưởng của phương pháp này là:    ‐ xấp xỉ hàm đối  tượng  f(x) bằng một hàm bậc 2 p2(x) qua 3 điểm cho  trước    ‐ cập nhật 3 điểm này bằng cách  thay một  trong 3 điểm bằng cực  tiểu  của hàm p2(x)    Qua 3 điểm:   [ ] [ ] [ ]{ }0 0 1 1 2 2(x ,f(x ) , (x ,f(x ) , (x ,f(x )   x0 < x1 < x2  ta tìm được đa thức nội suy p2(x) và điểm có đạo hàm bằng zero:    − + − + −= = ⎡ ⎤− + − + −⎣ ⎦ 2 2 2 2 2 2 0 1 2 1 2 0 2 0 1 3 0 1 2 1 2 0 2 0 1 f (x x ) f (x x ) f (x x )x x 2 f (x x ) f (x x ) f (x x )         (1)  Đặc biệt, nếu các  điểm tìm được trước đây phân bố đều với khoảng cách h(  nghĩa là x2 ‐ x1 = x1 ‐ x0 = h) thì (1) trở thành:    − + − + − − += = + − + −⎡ ⎤− + − + −⎣ ⎦ 2 2 2 2 2 2 0 1 2 1 2 0 2 0 1 0 1 2 3 0 0 1 20 1 2 1 2 0 2 0 1 f (x x ) f (x x ) f (x x ) 3f 4f fx x h 2( f 2f f )2 f (x x ) f (x x ) f (x x )   (2)  Ta  cập  nhật  3  điểm  theo  cách  này  cho  đến  khi  − ≈2 0x x 0   hay  − ≈2 0f(x ) f(x ) 0  và cực tiểu là x3. Quy tắc cập nhật 3 điểm là:    ‐ Trong trường hợp  < <0 3 1x x x  ta dùng  0 3 1(x , x , x )  hay  3 1 2(x , x , x )  làm 3  điểm mới tuỳ theo f(x3) < f(x1) hay không    ‐ Trong trường hợp  < <1 3 2x x x  ta dùng  1 3 2(x , x , x )  hay  0 1 3(x , x , x )  làm 3  điểm mới tuỳ theo f(x3) ≤ f(x1) hay không.  Quá trình tìm cực tiểu được mô tả trên hình sau:  375 Ta xây dựng hàm optquad() để thực hiện thuật toán này.  function [xo, fo] = optquad(f, x0, tolx, tolfun, maxiter)  %Tim cuc tieu cua f(x) bang phuong phap xap xi bac 2  if length(x0) > 2      x012 = x0(1:3);  else      if length(x0) == 2          a = x0(1);           b = x0(2);      else           a = x0 ‐ 10;           b = x0 + 10;      end      x012 = [a (a + b)/2  b];  end  f012 = f(x012);  [xo, fo] = optquad0(f, x012, f012, tolx, tolfun, maxiter);  function [xo, fo] = optquad0(f, x012, f012, tolx, tolfun, k)  x0 = x012(1);   x1 = x012(2);   x2 = x012(3);  376 f0 = f012(1);   f1 = f012(2);   f2 = f012(3);  nd = [f0 ‐ f2 f1 ‐ f0 f2 ‐ f1]*[x1*x1 x2*x2 x0*x0; x1 x2 x0]ʹ;  x3 = nd(1)/2/nd(2);   f3 = feval(f, x3); %Pt.(1)  if k <= 0 | abs(x3 ‐ x1) < tolx | abs(f3 ‐ f1) < tolfun      xo = x3;       fo = f3;      if k == 0          fprintf(ʹDay co the xem la diem cuc tieuʹ)      end  else      if x3 < x1          if f3 < f1               x012 = [x0 x3 x1];               f012 = [f0 f3 f1];          else               x012 = [x3 x1 x2];               f012 = [f3 f1 f2];          end      else          if f3 <= f1              x012 = [x1 x3 x2];               f012 = [f1 f3 f2];          else               x012 = [x0 x1 x3];               f012 = [f0 f1 f3];          end      end      [xo, fo] = optquad0(f, x012, f012, tolx, tolfun, k ‐ 1);  end  Để tìm điểm cực tiểu ta dùng chương trình ctoptquad.m:  clear all, clc  377 %f = inline(ʹ(x.^2 ‐ 4).^2/8 ‐ 1ʹ);  a = 0;   b = 3;  delta = 1e‐8;  epsilon = 1e‐10;  maxiter = 100;  [xmin, ymin] = optquad(f, [a  b], delta, epsilon, maxiter)  §4. PHƯƠNG PHÁP NELDER ‐ MEAD    Phương pháp Nelder ‐ Mead có thể dùng để tìm cực tiểu của hàm nhiều  biến mà phương pháp  tiết diện vàng hay phương pháp xấp xỉ bậc 2 không  dùng được. Thuật toán Nelder ‐ Mead gồm các bước như sau:  Bước 1: Cho 3 điểm đầu tiên a, b, c với f(a) < f(b) < f(c)  Bước 2: Nếu 3 điểm và các giá trị tương ứng của hàm đủ gần nhau thì ta  coi a là điểm cực tiểu và kết thúc quá trình tính  Bước  3:  Nếu  không  ta  coi  điểm  cực  tiểu  nằm  đối  diện  với  điểm  c  trên  đường  ab(xem  hình  vẽ) và lấy:    e = m + 2(m ‐ c)   với   m = (a + b)/2  và nếu f(e) < f(b) thì lấy:    r = (m + e)/2 = 2m ‐ c  và nếu f(r) <  f(c)  thì  lấy r  làm giá  trị mới  của  c;  nếu  f(r)  ≥  f(b)  thì  lấy:    s = (c + m)/2  và nếu f(s) < f(c) thì lấy s làm giá trị mới của c; nếu không bỏ các điểm b, c và  dùng m và c1 = (a + c)/2 làm điểm b và c mới và cho rằng cực tiểu nằm quanh  điểm a.    Bước 4: Trở về bước 1  Ta xây dựng hàm neldermead() để thực hiện thuật toán này:  function [xo, fo] = neldermead(f, x0, tolx, tolfun, maxiter)  n = length(x0);  e c1  a m r  c b s2  s1  c2  m = (a + b)/2  r = m + (m ‐ c)  e = m + 2(m ‐ c)  s1 = (c + m)/2  s2 = (m + r)/2  c1 = (c + a)/2  c2 = (r + a)/2  378 if n == 1 %truong hop ham 1 bien      [xo,fo] = optquad(f, x0, tolx, tolfun);       return  end  S = eye(n);  for i = 1:n       i1 = i + 1;       if i1 > n          i1 = 1;       end      abc = [x0; x0 + S(i,:); x0 + S(i1,:)];       fabc = [feval(f, abc(1, :)); feval(f,abc(2, :)); feval(f,abc(3, :))];      [x0, fo] = neldermead0(f, abc, fabc, tolx, tolfun, maxiter);      if n < 3,           break;       end   end  xo = x0;  function [xo, fo] = neldermead0(f, abc, fabc, tolx, tolfun, k)  [fabc, I] = sort(fabc);   a = abc(I(1), :);   b = abc(I(2), :);   c = abc(I(3), :);  fa = fabc(1);   fb = fabc(2);   fc = fabc(3);   fba = fb ‐ fa;   fcb = fc ‐ fb;  if k <= 0 | abs(fba) + abs(fcb) < tolfun | abs(b ‐ a) + abs(c ‐ b) < tolx      xo = a;       fo = fa;      if k == 0          fprintf(ʹXem day la diem cuc tieuʹ)      end  else  379     m = (a + b)/2;       e = 3*m ‐ 2*c;       fe = feval(f, e);      if fe < fb          c = e;           fc = fe;      else          r = (m + e)/2;           fr = feval(f, r);          if fr < fc              c = r;               fc = fr;           end          if fr >= fb              s = (c + m)/2;               fs = feval(f, s);              if fs < fc                  c = s;                   fc = fs;              else                   b = m;                   c = (a + c)/2;                   fb = feval(f, b);                   fc = feval(f, c);              end          end      end      [xo, fo] = neldermead0(f, [a;b;c], [fa fb fc], tolx, tolfun, k ‐ 1);  end  Để  tìm  cực  tiểu  của  hàm  = = − − + −2 21 1 2 1 2 2z f(x,y) x x x 4x x x   lân  cận  [0  0]  ta  dùng chương trình ctneldermead.m:  clear all, clc  f = inline(ʹx(1)*x(1) ‐ x(1)*x(2) ‐ 4*x(1) + x(2) *x(2) ‐ x(2)ʹ);  x0 = [0  0];   380 b = 1;  delta = 1e‐8;  epsilon = 1e‐10;  maxiter = 100;  [xmin, ymin] = neldermead(f, x0, delta, epsilon, maxiter)  §5. PHƯƠNG PHÁP ĐỘ DỐC LỚN NHẤT    Phương  pháp  này  tìm  điểm  cực  tiểu  của  hàm  n  biến  theo  hướng  gradient âm:    [ ] [ ] ⎡ ⎤∂ ∂ ∂− = −∇ = −⎢ ⎥∂ ∂ ∂⎣ ⎦L1 2 n f(x) f(x) f(x)g( x ) f( x ) x x x với kích thước bước tính αk tại lần lặp thứ k được hiệu chỉnh sao cho giá trị  hàm cực tiểu theo hướng tìm. Thuật toán gồm các bước sau:    • Tại lần lặp thứ k = 0, tìm giá trị hàm f(x0) với điểm khởi đầu x0  • Tại lần lặp thứ k, tìm αk theo hướng ‐g(x)    ( )− − − −α = − αk 1 k 1 k 1 k 1f x g / g             (1)  • Tính giá trị xk:    − − − −= − αk k 1 k 1 k 1 k 1x x g / g             (2)  • Nếu xk ≈ xk‐1 và f(xk) ≈ f(xk‐1) thì coi là cực tiểu, nếu không thì quay về  bước 2.    function [xo, fo] = steepest(f, x0, tolx, tolfun, alpha0,maxiter)  if nargin < 6      maxiter = 100;   end   if nargin < 5      alpha0 = 10; %kich thuoc khoi gan  end  if nargin < 4      tolfun = 1e‐8;   end %|f(x)| < tolfun mong muon  if nargin < 3      tolx = 1e‐6;   end %|x(k)‐ x(k ‐ 1)| < tolx mong muon  x = x0;   381 fx0 = feval(f, x0);   fx = fx0;  alpha = alpha0;   kmax1 = 25;  warning = 0;   for k = 1:maxiter      g = grad(f, x);       g = g/norm(g); %gradient la vec to hang      alpha = alpha*2; %de thu di theo huong gradient am      fx1 = feval(f, x ‐ alpha*2*g);      for k1 = 1:kmax1 %tim buoc toi uu          fx2 = fx1;        fx1 = feval(f, x ‐ alpha*g);          if fx0 > fx1+ tolfun & fx1  fx1 < fx2              den = 4*fx1 ‐ 2*fx0 ‐ 2*fx2;     num = den ‐ fx0 + fx2; %              alpha = alpha*num/den;Pt.(1)              x = x ‐ alpha*g;               fx = feval(f,x); %Pt.(2)              break;          else               alpha = alpha/2;          end      end      if k1 >= kmax1          warning = warning + 1; %kg tim duoc buoc toi uu      else           warning = 0;      end      if warning >= 2|(norm(x ‐ x0) < tolx&abs(fx ‐ fx0) < tolfun)          break;       end      x0 = x;       fx0 = fx;  end  xo = x; fo = fx;  382 if k ==maxiter      fprintf(ʹDay la ket qua tot nhat sau  %d lan lapʹ, maxiter)  end  function g = grad(func, x)  % Tinh gradient cua ham f(x).  h = 1.0e‐6;  n = length(x);  g = zeros(1, n);  f0 = feval(func, x);  for i = 1:n      temp = x(i);      x(i) = temp + h;      f1 = feval(func, x);      x(i) = temp;      g(1, i) = (f1 ‐ f0)/h;  end  Để tìm cực tiểu của hàm ta dùng chương trình ctsteepest.m:  clear all, clc  f = inline(ʹx(1)*x(1) ‐ x(1)*x(2) ‐ 4*x(1) + x(2) *x(2) ‐ x(2)ʹ);  x0 = [0.5  0.5];  tolx = 1e‐4;   tolfun = 1e‐9;   alpha0 = 1;   maxiter = 100;  [xo, fo] = steepest(f, x0, tolx, tolfun, alpha0, maxiter)  §6. PHƯƠNG PHÁP NEWTON    Việc tìm điểm cực tiểu của hàm f(x) tương đương với việc xác định x để  cho gradient g(x) của hàm f(x) bằng zero. Nghiệm của g(x) = 0 có thể tìm được  bằng cách dùng phương pháp Newton cho hệ phương trình phi tuyến. Hàm  newtons(x) dùng để tìm nghiệm của phương trình g(x) = 0 là:  function [x, fx, xx] = newtons(f, x0, tolx, maxiter)  383 h = 1e‐4;   tolfun = eps;   EPS = 1e‐6;  fx = feval(f, x0);  nf = length(fx);   nx = length(x0);  if nf ~= nx      error(ʹKich thuoc cua g va x0 khong tuong thich!ʹ);   end  if nargin < 4      maxiter = 100;   end  if nargin < 3      tolx = EPS;   end  xx(1, :) = x0(:).ʹ;   for k = 1: maxiter      dx = ‐jacob(f, xx(k, :), h)\fx(:);%‐[dfdx]ˆ‐1*fx      xx(k + 1, :) = xx(k, :) + dx.ʹ;      fx = feval(f, xx(k + 1, :));       fxn = norm(fx);      if fxn < tolfun | norm(dx) < tolx          break;       end  end  x = xx(k + 1, :);  if k == maxiter      fprintf(ʹKet qua tot nhat sau %dlan lap\nʹ, maxiter)  end  function g = jacob(f, x, h) %Jacobian cua f(x)  if nargin < 3      h = 1e‐4;   end  h2 = 2*h;   n = length(x);   384 x = x(:).ʹ;   I = eye(n);  for n = 1:n      g(:, n) = (feval(f, x + I(n, :)*h) ‐ feval(f, x ‐ I(n,:)*h))ʹ/h2;  end  Để tìm cực tiểu của hàm bằng phương pháp Newtons ta dùng chương trình  ctnewtons.m:  clear all, clc  f = inline(ʹx(1).^2 ‐ x(1)*x(2) ‐ 4*x(1) + x(2).^2 ‐ x(2)ʹ);  g = inline(ʹ[(2*x(1) ‐ x(2) ‐4)  ( 2*x(2) ‐ x(1) ‐ 1)]ʹ);  x0 = [0.1  0.1];  tolx = 1e‐4;   maxiter = 100;  [xo, fo] = newtons(g, x0, tolx, maxiter)  §7. PHƯƠNG PHÁP GRADIENT LIÊN HỢP  1. Khái niệm chung: Một trong các phương pháp giải bài tón tìm cực tiểu của  hàm nhiều biến  là tìm cực  tiểu theo một biến  liên tiếp để đến gần điểm cực  tiểu. Chiến thuật chung là:    • chọn điểm [x0]  • lặp với i = 1, 2, 3,...    ‐ chọn vec tơ [vi]  ‐ cực tiểu hoá f([x]) dọc theo đường [xi‐1] theo hướng [vi]. Coi [xi]   là cực tiểu. Nếu  [ ] [ ]−− < εi i 1x x  thì kết thúc lặp  • kết thúc   2. Gradient liên hợp: Ta khảo sát hàm bậc 2:    [ ]( ) [ ] [ ] [ ] [ ][ ]= − + = − +∑ ∑∑ T Ti i i ,j i j i i j 1 1f x c b x A x x c b x x A x 2 2   (1)  đạo hàm của hàm theo xi cho ta:    ∂ = − +∂ ∑i i ,j jji f b A x x Viết dưới dạng vec tơ ta có:    [ ] [ ][ ]∇ = − +f b A x                 (2)  385 với ∇f là gradient của f.    Bây giờ  ta khảo  sát  sự  thay đổi gradient khi  ta di chuyển  từ  [x0]  theo  hướng của vec tơ [u] dọc theo đường:    [x] = [x0] + s[u]  với s là khoảng cách di chuyển. Thay biểu thức này vào (2) ta có gradient dọc  theo [u]:      [ ] [ ] [ ] [ ] [ ] [ ]( ) [ ] [ ][ ]+∇ = − + + = ∇ +0 00x s u xf b A x s u f s A u   Như vậy sự thay đổi gradient là s[A][u]. Nếu ta di chuyển theo hướng vuông  góc với vec tơ [v], nghĩa là   [v]T[u] = 0 hay [v]T[A][u] = 0              (3)  thì hướng của [u] và [v] là liên hợp. Điều này cho thấy khi ta muốn cực tiểu  hoá f(x) theo hướng [v], ta cần di chuyển theo hướng [u] để không làm hỏng  cực tiểu trước đó. Với hàm bậc hai n biến độc lập ta có thể xây dựng n hướng  liên hợp. Các phương pháp gradient  liên hợp khác nhau dùng các kỹ  thuật  khác nhau để tìm hướng liên hợp.  3. Phương pháp Powell: Phương pháp Powell  là phương pháp bậc zero, chỉ  đòi hỏi tính f([x]). Thuật toán gồm các bước:    • chọn điểm [x0]  •  chọn vec  tơ  [vi],  thường  lấy  [vi]  =  [ei] với  [ei]  là vec  tơ  đơn vị  theo  hướng [xi]  • vòng tròn    ‐ lặp với i = 1, 2,...  ∗  tìm  cực  tiểu  của  f([x]) dọc  theo  đường  đi  qua  [xi‐1]  theo  hướng [vi]; coi [xi] là cực tiểu       ‐ kết thúc lặp  ‐ [vn+1] = [x0] ‐ [xn]; tìm cực tiểu của f([x]) dọc theo đường đi qua[x0]  theo hướng [vn+1]; coi [xn+1] là cực tiểu  ‐ nếu  [ ] [ ]+ − < εn 1 nx x  thì kết thúc vòng lặp    ‐ lặp  ∗ [vi+1] = [v]    • kết thúc vòng tròn  Ta xây dựng hàm powell() để thực hiện thuật toán trên:  function [xmin, fmin, ncyc] = powell(h, tol)  % Phuong phap Powell de tim cuc tieu cua ham f(x1,x2,...,xn).  386 global x func v  if nargin < 2;       tol = 1.0e‐6;   end  if nargin < 1       h = 0.1;   end  if size(x,2) > 1      x = xʹ;   end   n = length(x);   df = zeros(n, 1);   u = eye(n);       xold = x;      fold = feval(func, xold);      for i = 1:n          v = u(1:n, i);          [a, b] = goldbracket(@fline, 0.0, h);          [s, fmin] = golden(@fline, a, b);          df(i) = fold ‐ fmin;          fold = fmin;          x = x + s*v;      end      v = x ‐ xold;  [a, b] = goldbracket(@fline, 0.0, h);      [s,fMin] = golden(@fline, a, b);      x = x + s*v;      if sqrt(dot(x‐xold, x‐xold)/n) < tol          xmin = x;           ncyc = j;           return      end      imax = 1;       dfmax = df(1);      for i = 2:n          if df(i) > dfmax  387             imax = i;               dfmax = df(i);          end      end      for i = imax:n‐1          u(1:n, i) = u(1:n, i+1);      end      u(1:n, n) = v;      end  error(ʹPhuong phap Powell khong hoi tuʹ)  function z = fiine(s) % f theo huong cua v  global x func v  z = feval(func, x+s*v);    function y = f(x)  y = 100.0*(x(2) ‐ x(1).^2).^2 + (1.0 ‐ x(1)).^2;  Để tìm điểm cực tiểu ta dùng chương trình ctpowell.m:  clear all, clc  global x  func   func = @f;  x = [‐1.0; 1.0];  [xmin, fmin, numcycles] = powell  fminsearch(func, x)  3. Phương pháp Fletcher ‐ Reeves: Phương pháp Powell cần n đường cực tiểu  hoá. Ta có thể chỉ cần dùng một đường với phương pháp bậc 1. Phương pháp  này có 2 phương án:  thuật  toán Fletcher  ‐ Reeves(FR) và  thuật  toán Polak  ‐  Ribiere(PR). Thuật toán tóm tắt như sau:    ‐ cho [x0], tính f([x0])     ‐ khở