Минимизация функции градиентным методом

[alenov]

Добрый день,
Нужно минимизировать нелинейную функцию при ряде ограничений. На данный момент данную задачу удалось решить в Excel с помощью модуля Solver с использованием метода сопряженных градиентов. Нужно воспроизвести это решение в моей программе. Я включил в целевую функцию ограничивающие условия, помноженные на множители Лагранжа. Таким образом, решается задача безусловной минимизации. Начиная из определенной точки, рассчитывается градиент и подбирается такая длина шага, чтобы он был максимальным и при этом значение целевой функции было минимальным, меньше, чем на предыдущем. Проблема в том, что в Excel значение целевой функции не обязательно уменьшается на каждом из шагов, а может и увеличиваться. Я не понимаю, каким образом это может происходить – в известных мне описаниях метода сопряженных градиентов значение функции на каждом последующем шаге должно быть меньше, чем на предыдущем. Не может ли кто-нибудь помочь разобраться как работает алгоритм в Excel?

[мат-ламер]

alenov в сообщении #776555 писал(а):

Проблема в том, что в Excel значение целевой функции не обязательно уменьшается на каждом из шагов, а может и увеличиваться.

Непонятно, почему это для Вас проблема? Вы свою программу пишите? Или хотите полную копию алгоритма Excela?

alenov в сообщении #776555 писал(а):

Я включил в целевую функцию ограничивающие условия, помноженные на множители Лагранжа. Таким образом, решается задача безусловной минимизации.

Любопытно. Откуда такой подход? Обычно штрафную функцию тоже добавляют.

[alenov]

Думаю, нужно больше деталей.
Нужно минимизировать функцию Гиббса-Дюхема:
$\min G(x)=\sum(x_i(G_0_i+RT\ln(x_i/\sum(x_i))))$, $G_0_i = \operatorname{const}, R = 1.985877534 cal, T = 298.15°K$
при $\sum(a_i_j x_i)=b_j; x_i>0$

$x_i$ - ?

После включения условий, помноженных на множители Лагранжа, и замены $x_i$ на $\exp(x'_i)$, чтобы выполнить условие положительности переменных, целевая функция принимает вид:

$\min G(x)=\sum(\exp(x'_i)(G_0_i /RT+x'_i-\ln(\sum(\exp(x'_i)))))+\sum(L_j b_j(\sum(a_i_j \exp(x'_i))))$

на шаге k значение переменной рассчитывается из $x'^{(k)}_i -step \cdot dG^{(k)}/dx'_i$, т.е. задача сводится к одномерной минимизации целевой функции по переменной t, в настоящее время я пытаюсь применить метод золотого сечения.

Производные рассчитываются по формулам:
$dG/dx'_i=G_0_i /RT+x'_i-\ln(\sum(\exp(x'_i)))+\sum(L_j a_i_j)$ и
$dG/dL_j=\sum(a_i_j \exp(x'_i)-b_j)$

[ИСН]

Химические функции обычно просты и выпуклы, без ложных минимумов. Так что не надо умных методов. Умные заведут в болото. Тут должны работать тупые.

[alenov]

Копия алгоритма используемого в Excel не нужна, нужно просто найти минимум этой функции, желательно разными методами, чтобы выбрать оптимальный. Работаю над своей программой. В принципе не важно как решить, главное решить. Excel находит минимум очень хорошо, вот и пытаюсь изучить, как он это делает. Можно и штрафными функциями. Но как мне кажется, они больше для ограничений выраженных в неравенствах подходят.

-- 17.10.2013, 21:12 --

Что вы подразумеваете под тупыми методами? Симплекс-метод? Он вроде не подходит, т.к. есть логарифм - целевая функция нелинейная.

-- 17.10.2013, 21:24 --

В болоте я уже нахожусь Какие максимально тупые методы вы можете посоветовать? Чем проще, тем лучше.

[мат-ламер]

alenov в сообщении #776581 писал(а):

В болоте я уже нахожусь Какие максимально тупые методы вы можете посоветовать? Чем проще, тем лучше.

Если Вы в болоте, то лучше воспользоваться готовыми программами. Если Excel не подходит, то воспользуетесть МатЛабом.

[alenov]

К сожалению, excel и matlab не подходят, т.к. этот расчет часть большого алгоритма. Нет возможности проводить эту операцию в сторонних программах. Если вы посоветуете open source программу для нахождения минимума этой функции, буду очень признателен. Понятно, что это сложно, но нет другого выхода, надо разобраться с алгоритмом, уже не первый месяц этим занимаюсь. В данный момент у меня есть конкретный вопрос: что минимизируется в excel, если целевая функция может возрастать?

Пока тема закрыта, напишу здесь.
Продолжаю работать над своей задачей.
Очень благодарен Мат-Ламеру за указание использовать штрафные функции. Раньше я думал, что они мне не подходят. С использованием штрафных функций у меня кое-что получилось.

Сконструировал вот такую целевую функцию:
$ \min P(x^{(k)}, r^{(k)})=\sum(x^{(k)}_i(G_0_i+RT\ln(x^{(k)}_i/\sum(x^{(k)}_i))))+0.5r^{(k)} \sum(h_j(x^{(k)})^2)-r^{(k)} \sum(\ln(-x^{(k)}_i)) $ (последний терм не равен 0 только, если $ x^{(k)}_i<0$)
при $ h_j(x^{(k)})=\sum (b_j-a_j_ix^{(k)}_i )$

Минимизирую ее на шаге k+1 с помощью метода покоординатного спуска (самый простой метод): $ x^{(k+1)}_i=x^{(k)}_i-step\cdot g^{(k)}_i / |g^{(k)}_i| $, где $g^{(k)}_i=dP(x^{(k)})/dx_i$. Она сходится к минимуму, рассчитанному в Excel, но не достаточно точно, процентов на 10 отличается.
При поиске минимума у меня возникла следующая проблема. Значения функций Гиббса-Дюхема и штрафной функции существенно различается, Гиббса-Дюхема на 3-4 порядка больше. Таким образом, частные производные не могут изменить свой знак практически никогда. Я думаю, что функции нужно масштабировать, но как это правильно сделать, пока не могу найти. Пока я просто делю Гиббса-Дюхема на 1000, но это конечно неверно, хоть и дает некий прогресс в решении задачи. Значение r у меня пока всегда единица. Вопрос: как правильно масштабировать функции?

[Deggial]

Тема перемещена из форума «Помогите решить / разобраться (М)» в форум «Карантин»
Причина переноса: формулы плохо оформлены

alenov, наберите все формулы и термы $\TeX$ом нормально. Инструкции по оформлению формул здесь или здесь (или в этом видеоролике). Программный код можете оформлять тегом code.
После исправлений сообщите в теме Сообщение в карантине исправлено, и тогда тема будет возвращена.

Тема перемещена из форума «Карантин» в форум «Помогите решить / разобраться (М)»
вернул

[alenov]

Deggial, большое спасибо!

На данный момент главный вопрос: как правильно масштабировать целевую функцию и/или ее производные?

[мат-ламер]

alenov в сообщении #776599 писал(а):

Очень благодарен Мат-Ламеру за указание использовать штрафные функции.

Я не предлагал использовать штрафные функции. Я говорил, что обычно к функции Лагранжа добавляют штрафные функции. Получается метод модифицированной функции Лагранжа. Там возможно и масштабировать не надо. Градиентный метод практически неработоспособен. МатЛаб можно стыковать с внешними программами. И наоборот.

alenov в сообщении #776599 писал(а):

В данный момент у меня есть конкретный вопрос: что минимизируется в excel, если целевая функция может возрастать?

Что понимать под целевой функцией? Когда мы задачу с ограничениями сводим к задаче без ограничений, целевая функция меняется. Новая будет уменьшаться. А старая может и возрастать.

[alenov]

После большого перерыва хочу опять попросить вашей помощи. Указания Мат-Ламера мне очень помогли, я сделал как он сказал - воспользовался методом комбинирования минимизируемой функции с множителями Лагранжа и с штрафной функцией. Я написал код, который решает задачу. Но есть проблема. Искомое решение с заданной точностью ищется слишком медленно. На моем компьютере поиск занимает 35 секунд, в то время как требуется решать задачу за 0.1 секунды. Пока код не оптимизирован, но дело не в этом, а в математике (пробовал заняться оптимизацией, но выйти из 30 секунд все равно не получается). За первые 2 цикла приближение идет очень быстро (точность 1е-4 (значение штрафной функции)). А вот потом код работает практически в холостую и требуемая точность (1е-10) достигается примерно за 1000 циклов. Для минимизации целевой функции я использую метод сопряженных градиентов Fletcher and Reeves (Bazaraa, 2006, страница 423). У меня не получилось сделать, как там написано, выход из цикла если градиент целевой функции меньше критического значения. У меня градиент не уменьшается, а колеблется (то больше, то меньше, но почти никогда не бывает очень маленьким). Косвенное подтверждение этому феномену я нашел в этой книге (Бертсекас, 1982, страница 66). Но может быть я просто что-то делаю неправильно или можно применить какой-то другой более эффективный метод для минимизации по направлению. В моем алгоритме выход происходит, когда целевая функция перестает существенно изменяться. Я использую довольно странный метод одномерной минимизации, который я взял из (Capitani and Brown, 1987), но для этой задачи он в несколько раз эффективнее метода золотого сечения. В общем можно сказать, что я много чего перепробовал и уже даже и не знаю, что можно еще сделать, чтобы повысить эффективность вычислений.
Вот собственно код, написанный на Delphi, выход из которого нужно сделать более эффективным:

repeat
    dk_max:=-1e+100;
    for i := 0 to n-1 do
     if abs(dk[i])>dk_max then dk_max:=abs(dk[i]);
    step:=1/dk_max;
    Gbulk_old:=Gbulk;
    k:=0;
    repeat                  //One-dimentional minimization (Capitani and Brown, 1987)
     for i := 0 to n-1 do
      x_1[i]:=x[i]+step*dk[i];
     Gbulk_1:=GibbsDuhem(x_1,G0,n)+Lagrange(x_1,Lag,aij,b,n,m)+Penalty(x_1,aij,b,rk,n,m)+Lagrange_pH(x_1,Lag,aij,b,n,m)+Penalty_pH(x_1,aij,b,rk,n,m);
     if Gbulk_1<Gbulk then begin
      step:=step+step;
      for i := 0 to n-1 do
       x[i]:=x_1[i];
      Gbulk:=Gbulk_1;
     end else begin
      f:=1;
      repeat
       step:=step/2/f;
       for i := 0 to n-1 do
        x_1[i]:=x[i]+step*dk[i];
       Gbulk_1:=GibbsDuhem(x_1,G0,n)+Lagrange(x_1,Lag,aij,b,n,m)+Penalty(x_1,aij,b,rk,n,m)+Lagrange_pH(x_1,Lag,aij,b,n,m)+Penalty_pH(x_1,aij,b,rk,n,m);
       Inc(f);
      until (Gbulk_1<Gbulk) or (f>100) or (abs(step*dk_max)<eps);
     end;
     Inc(k);
    until (Gbulk_1>Gbulk) or (k>100);       //Finish of one-dimentional minimization

    if Gbulk_old-Gbulk<eps*abs(Gbulk) then begin    //Exit from the cycle if the function doesn't change
     if lll=0 then
      for i := 0 to n-1 do
       x_bu[i]:=x[i];
     Inc(lll);
     if lll>10 then break;
    end else lll:=0;

    gradx_old:=gradx;
    gradx:=0;
    for i:=0 to n-1 do begin                        //Calculate new gradient
     g[i]:=GibbsDuhemDerivative(x,G0,i,n)+LagrangeDerivative(Lag,aij,i,m)+PenaltyDerivative(x,aij,b,rk,i,n,m)+PenaltyDerivative_pH(x,aij,b,rk,i,n,m)+LagrangeDerivative_pH(Lag,aij,i,m);
     gradx:=gradx+g[i]*g[i];
    end;
    Betta:= gradx/gradx_old;                        //Calculate new conjugate directions
    if fff>20 then begin                            //Refresh of conjugate directions every 20 steps
     for i:=0 to n-1 do
      dk[i]:=-g[i];
     fff:=0; end else begin
     for i:=0 to n-1 do
      dk[i]:=-g[i]+Betta*dk[i];
     Inc(fff);
    end;
    Inc(o);
   until o>1000;

Я не знаю, могу ли я вставлять сюда такие большие куски кода, но возможно, что неэффективная процедура зашита в другом месте. Поэтому рискну вставить сюда весь текст програмы. Также исходники можно скачать здесь (сама програма в папке /Win32/Debug, она берет исходные параметры для вычислений из файла ini): https://drive.google.com/file/d/0B08REaASOrMBSVd3RDRqWk1hOUk/edit?usp=sharing

unit Lagrange_Multipliers;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics,
  Controls, Forms, Dialogs, Grids, StdCtrls, IniFiles, Math, DateUtils;

type
  TA = array of array of extended;
  TForm1 = class(TForm)
    StringGrid1: TStringGrid;
    StringGrid2: TStringGrid;
    EPSEdit: TEdit;
    TimeEdit: TEdit;
    Label1: TLabel;
    Label2: TLabel;
    Button14: TButton;
    procedure FormClose(Sender: TObject; var Action: TCloseAction);
    procedure FormPaint(Sender: TObject);
    procedure Button14Click(Sender: TObject);

  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;
  IniFileName: string;
  IniFile: TIniFile;
  IniFormSaveName: string;
  IniFormSaveFile: TIniFile;

const
  R: extended = 1.985877534;
  T: extended = 298.15;

implementation

{$R *.dfm}

procedure ReadIni;       //Read data from ini
var
 IniPath: string;
 k, m: integer;
begin
 GetDir(0,IniPath);
 IniFileName:=IniPath+'\Crono.ini';
  if FileExists(IniFileName) then
   begin
    IniFile:=TIniFile.Create(IniFileName);
    Form1.Left:=IniFile.ReadInteger('position','left',100);
    Form1.Top:=IniFile.ReadInteger('position','top',100);
    Form1.StringGrid1.ColWidths[0]:=20;
    Form1.StringGrid1.Cells[1,0]:='Phase';
    Form1.StringGrid1.ColWidths[1]:=60;
    Form1.StringGrid1.Cells[2,0]:='G';
    Form1.StringGrid1.ColWidths[2]:=60;
    Form1.StringGrid1.Cells[3,0]:='O';
    Form1.StringGrid1.ColWidths[3]:=20;
    Form1.StringGrid1.Cells[4,0]:='H';
    Form1.StringGrid1.ColWidths[4]:=20;
    Form1.StringGrid1.Cells[5,0]:='Si';
    Form1.StringGrid1.ColWidths[5]:=20;
    Form1.StringGrid1.Cells[6,0]:='z';
    Form1.StringGrid1.ColWidths[6]:=20;
    Form1.StringGrid1.Cells[7,0]:='pH';
    Form1.StringGrid1.ColWidths[7]:=20;
    Form1.StringGrid1.Cells[8,0]:='x';
    Form1.StringGrid2.Cells[0,0]:='Phase';
    Form1.StringGrid2.ColWidths[0]:=35;
    Form1.StringGrid2.Cells[0,1]:='O';
    Form1.StringGrid2.Cells[0,2]:='H';
    Form1.StringGrid2.Cells[0,3]:='Si';
    Form1.StringGrid2.Cells[0,4]:='z';
    Form1.StringGrid2.Cells[0,5]:='pH';
    Form1.StringGrid2.Cells[0,6]:='G';
    Form1.StringGrid2.Cells[0,7]:='P';
    Form1.StringGrid2.Cells[1,0]:='a*x';
    Form1.StringGrid2.Cells[2,0]:='g';
    Form1.StringGrid2.Cells[3,0]:='Lag';
    Form1.StringGrid2.Cells[4,0]:='InSum';
    Form1.StringGrid2.ColWidths[4]:=60;
    m:=1;
    for k := 1 to Form1.StringGrid1.RowCount-1 do
     begin
      Form1.StringGrid1.Cells[0,k]:=IntToStr(m);
      Form1.StringGrid1.Cells[1,k]:=IniFile.ReadString('Name of phase',IntToStr(m),'');
      Form1.StringGrid1.Cells[2,k]:=IniFile.ReadString('G',IntToStr(m),'');
      Form1.StringGrid1.Cells[3,k]:=IniFile.ReadString('O',IntToStr(m),'');
      Form1.StringGrid1.Cells[4,k]:=IniFile.ReadString('H',IntToStr(m),'');
      Form1.StringGrid1.Cells[5,k]:=IniFile.ReadString('Si',IntToStr(m),'');
      Form1.StringGrid1.Cells[6,k]:=IniFile.ReadString('z',IntToStr(m),'');
      Form1.StringGrid1.Cells[7,k]:=IniFile.ReadString('pH',IntToStr(m),'');
      Inc(m);
     end;
    m:=1;
    for k := 1 to Form1.StringGrid2.RowCount-1 do
     begin
      Form1.StringGrid2.Cells[4,k]:=IniFile.ReadString('InSum',IntToStr(m),'');
      Inc(m);
     end;
    IniFile.Free;
   end else
  Form1.Position:=poScreenCenter;
end;

procedure WriteIni;       //Save data in ini
var
 k, m: integer;
begin
 if FileExists(IniFileName) then
  if not(DeleteFile(IniFileName)) then
    ShowMessage('Crono did not delete ini file');
 IniFile:=TIniFile.Create(IniFileName);
 IniFile.WriteInteger('position','left',Form1.Left);
 IniFile.WriteInteger('position','top',Form1.Top);
 m:=1;
 for k := 1 to Form1.StringGrid1.RowCount-1 do
  begin
   IniFile.WriteString('Name of phase',IntToStr(m),Form1.StringGrid1.Cells[1,k]);
   IniFile.WriteString('G',IntToStr(m),Form1.StringGrid1.Cells[2,k]);
   IniFile.WriteString('O',IntToStr(m),Form1.StringGrid1.Cells[3,k]);
   IniFile.WriteString('H',IntToStr(m),Form1.StringGrid1.Cells[4,k]);
   IniFile.WriteString('Si',IntToStr(m),Form1.StringGrid1.Cells[5,k]);
   IniFile.WriteString('z',IntToStr(m),Form1.StringGrid1.Cells[6,k]);
   IniFile.WriteString('pH',IntToStr(m),Form1.StringGrid1.Cells[7,k]);
   Inc(m);
  end;
 m:=1;
 for k := 1 to Form1.StringGrid2.RowCount-1 do
  begin
   IniFile.WriteString('InSum',IntToStr(m),Form1.StringGrid2.Cells[4,k]);
   Inc(m);
  end;
 IniFile.Free;
end;

function GibbsDuhem(x:array of extended; G0:array of extended; n: integer): extended;
var
 xsum, Gbulk: extended;
 i: integer;
begin
 xsum:=0;
 for i := 0 to n-1 do
  xsum:=xsum+exp(x[i]);
 Gbulk:=0;
 for i := 0 to n-1 do
  Gbulk:=Gbulk+exp(x[i])*(G0[i]+x[i]-ln(xsum));
 Result:=Gbulk;
end;

function Lagrange(x: array of extended; Lag: array of extended; aij: TA; b: array of extended; n: integer; m: integer): extended;
var
 P, ggg: extended;
 i, k: integer;
begin
 P:=0;
 for k := 0 to m-2 do begin
  ggg:=0;
  for i := 0 to n-1 do
   ggg:=ggg+exp(x[i])*aij[k,i];
  P:=P+Lag[k]*(ggg-b[k]);
 end;
 Result:=P;
end;

function Lagrange_pH(x: array of extended; Lag: array of extended; aij: TA; b: array of extended; n: integer; m: integer): extended;
var
 xmult: extended;
 l: integer;
begin
 xmult:=0;
 for l := 0 to n-1 do
  xmult:=xmult+aij[m-1,l]*x[l];
 Result:=Lag[m-1]*(xmult/ln(10)+b[m-1]);                    //ln(10)=2.30258509299405
end;

function Penalty(x:array of extended; aij: TA; b: array of extended; rk: extended; n: integer; m: integer): extended;
var
 P, ggg: extended;
 i, k: integer;
begin
 P:=0;
 for k := 0 to m-2 do begin
  ggg:=0;
  for i := 0 to n-1 do
   ggg:=ggg+exp(x[i])*aij[k,i];
  P:=P+(ggg-b[k])*(ggg-b[k]);
 end;
 Result:=0.5*rk*P;
end;

function Penalty_pH(x:array of extended; aij: TA; b: array of extended; rk: extended; n: integer; m: integer): extended;
var
 xmult: extended;
 l: integer;
begin
 xmult:=0;
 for l := 0 to n-1 do
  xmult:=xmult+aij[m-1,l]*x[l];
 Result:=0.5*rk*(xmult/ln(10)+b[m-1])*(xmult/ln(10)+b[m-1]);                  //ln(10)=2.30258509299405
end;

function GibbsDuhemDerivative(x:array of extended; G0:array of extended; i: integer; n:integer): extended;
var
 xsum: extended;
 l: integer;
begin
 xsum:=0;
 for l := 0 to n-1 do
  xsum:=xsum+exp(x[l]);
 Result:=G0[i]+x[i]-ln(xsum);
end;

function LagrangeDerivative(Lag: array of extended; aij: TA; i: integer; m:integer): extended;
var
 aaa: extended;
 l: integer;
begin
 aaa:=0;
 for l := 0 to m-2 do
  aaa:=aaa+aij[l,i]*Lag[l];
 Result:=aaa;
end;

function LagrangeDerivative_pH(Lag: array of extended; aij: TA; i: integer; m:integer): extended;
begin
 Result:=aij[m-1,i]*Lag[m-1]/ln(10);                 //ln(10)=2.30258509299405
end;

function PenaltyDerivative(x:array of extended; aij: TA; b: array of extended; rk: extended; i: integer; n:integer; m:integer): extended;
var
 aaa, ggg: extended;
 l, k: integer;
begin
 aaa:=0;
 for l := 0 to m-2 do begin
  ggg:=0;
  for k := 0 to n-1 do
   ggg:=ggg+exp(x[k])*aij[l,k];
  aaa:=aaa+aij[l,i]*(ggg-b[l]);
 end;
 Result:=rk*aaa;
end;

function PenaltyDerivative_pH(x:array of extended; aij: TA; b: array of extended; rk: extended; i: integer; n:integer; m:integer): extended;
var
 xmult: extended;
 l: integer;
begin
 xmult:=0;
 for l := 0 to n-1 do
  xmult:=xmult+aij[m-1,l]*x[l];
 Result:=aij[m-1,i]*rk/ln(10)*(xmult/ln(10)+b[m-1]);          //ln(10)=2.30258509299405
end;

procedure TForm1.Button14Click(Sender: TObject);
var
 i, k, l, m, n, o, ii, lll, f, fff: integer;
 eps, step, Gbulk, Gbulk_old, Gbulk_1, rk, gradx, gradx_old, Betta, P_new, P_old, ggg, dk_max: extended;
 aij: TA;
 G0, x, b, g, dk, x_1, x_bu, Lag: array of extended;
 TimeStart: TDateTime;
begin
 try
  TimeStart:=Now;
  n:=0;
  for i := 1 to Form1.StringGrid1.RowCount-1 do
   if Form1.StringGrid1.Cells[6,i]<>'' then Inc(n);
  SetLength(G0,n);
  m:=Form1.StringGrid1.ColCount-4;
  SetLength(aij,n+1,m+1);
  k:=0;
  for i := 1 to Form1.StringGrid1.RowCount-1 do
   if Form1.StringGrid1.Cells[6,i]<>'' then begin
    G0[k]:=StrToFloat(Form1.StringGrid1.Cells[2,i])/(R*T);       //Standard Gibbs energies
    Inc(k);
   end;
  l:=0;
  repeat
   k:=0;
   repeat
    aij[l,k]:=StrToFloat(StringGrid1.Cells[l+3,k+3]);     //Stoichiometric coefficients
    Inc(k);
   until k=n;
   Inc(l);
  until l=m;
  SetLength(x,n+1);                //Variables
  SetLength(x_1,n+1);
  SetLength(x_bu,n+1);
  for i:=0 to n-1 do
   x[i]:=0;
  SetLength(Lag, m+1);               //Lagrange Multipliers
  for i:=0 to m-1 do
   Lag[i]:=0;
  eps:=StrToFloat(EPSEdit.Text);       //Precision parameter
  rk:=100;                 //The initial rk value (100) was taken from (Lantagne et al., 1988)
  SetLength(b,m);                      //Constraints
  for i := 1 to m do
   b[i-1]:=StrToFloat(StringGrid2.Cells[4,i]);
  SetLength(g,n);                      //Gradient
  SetLength(dk,n);                     //Conjugate directions
  l:=0;
  repeat                        //Start of the minimization cycle
   o:=0;
   lll:=0;
   fff:=0;
   P_new:=Penalty(x,aij,b,rk,n,m)+Penalty_pH(x,aij,b,rk,n,m);
   Gbulk:=GibbsDuhem(x,G0,n)+Lagrange(x,Lag,aij,b,n,m)+Lagrange_pH(x,Lag,aij,b,n,m)+P_new;
   gradx:=0;
   for i:=0 to n-1 do begin
    g[i]:=GibbsDuhemDerivative(x,G0,i,n)+LagrangeDerivative(Lag,aij,i,m)+PenaltyDerivative(x,aij,b,rk,i,n,m)+PenaltyDerivative_pH(x,aij,b,rk,i,n,m)+LagrangeDerivative_pH(Lag,aij,i,m);
    dk[i]:=-g[i];
    gradx:=gradx+g[i]*g[i];
   end;
   repeat
    dk_max:=-1e+100;
    for i := 0 to n-1 do
     if abs(dk[i])>dk_max then dk_max:=abs(dk[i]);
    step:=1/dk_max;
    Gbulk_old:=Gbulk;
    k:=0;
    repeat                  //One-dimentional minimization (Capitani and Brown, 1987)
     for i := 0 to n-1 do
      x_1[i]:=x[i]+step*dk[i];
     Gbulk_1:=GibbsDuhem(x_1,G0,n)+Lagrange(x_1,Lag,aij,b,n,m)+Penalty(x_1,aij,b,rk,n,m)+Lagrange_pH(x_1,Lag,aij,b,n,m)+Penalty_pH(x_1,aij,b,rk,n,m);
     if Gbulk_1<Gbulk then begin
      step:=step+step;
      for i := 0 to n-1 do
       x[i]:=x_1[i];
      Gbulk:=Gbulk_1;
     end else begin
      f:=1;
      repeat
       step:=step/2/f;
       for i := 0 to n-1 do
        x_1[i]:=x[i]+step*dk[i];
       Gbulk_1:=GibbsDuhem(x_1,G0,n)+Lagrange(x_1,Lag,aij,b,n,m)+Penalty(x_1,aij,b,rk,n,m)+Lagrange_pH(x_1,Lag,aij,b,n,m)+Penalty_pH(x_1,aij,b,rk,n,m);
       Inc(f);
      until (Gbulk_1<Gbulk) or (f>100) or (abs(step*dk_max)<eps);
     end;
     Inc(k);
    until (Gbulk_1>Gbulk) or (k>100);       //Finish of one-dimentional minimization

    if Gbulk_old-Gbulk<eps*abs(Gbulk) then begin    //Exit from the cycle if the function doesn't change
     if lll=0 then
      for i := 0 to n-1 do
       x_bu[i]:=x[i];
     Inc(lll);
     if lll>10 then break;
    end else lll:=0;

    gradx_old:=gradx;
    gradx:=0;
    for i:=0 to n-1 do begin                        //Calculate new gradient
     g[i]:=GibbsDuhemDerivative(x,G0,i,n)+LagrangeDerivative(Lag,aij,i,m)+PenaltyDerivative(x,aij,b,rk,i,n,m)+PenaltyDerivative_pH(x,aij,b,rk,i,n,m)+LagrangeDerivative_pH(Lag,aij,i,m);
     gradx:=gradx+g[i]*g[i];
    end;
    Betta:= gradx/gradx_old;                        //Calculate new conjugate directions
    if fff>20 then begin                            //Refresh of conjugate directions every 20 steps
     for i:=0 to n-1 do
      dk[i]:=-g[i];
     fff:=0; end else begin
     for i:=0 to n-1 do
      dk[i]:=-g[i]+Betta*dk[i];
     Inc(fff);
    end;
    Inc(o);
   until o>1000;

   if lll>0 then
    for i := 0 to n-1 do
     x[i]:=x_bu[i];
   P_old:=P_new;
   P_new:=Penalty(x,aij,b,rk,n,m)+Penalty_pH(x,aij,b,rk,n,m);
   if (P_old>P_new) and (P_new>P_old/4) then
    rk:=rk*32;                  //The multiplier value (32) was taken from (Lantagne et al., 1988)
   for i := 0 to m-2 do begin               //Calculate Lagrange multipliers
    ggg:=0;
    for ii := 0 to n-1 do
     ggg:=ggg+exp(x[ii])*aij[i,ii];
    Lag[i]:=Lag[i]+rk*(ggg-b[i]);
   end;
   i:=m-1;
   ggg:=0;
   for ii := 0 to n-1 do
    ggg:=ggg+aij[i,ii]*x[ii];
   Lag[i]:=Lag[i]+rk*(ggg/ln(10)+b[i]);
   Inc(l);
  until (l>1000) or (sqrt(2*P_new/rk)<eps);           //Finish of the minimization cycle

  for i := 0 to n-1 do                                          //Write the obtained result in the table
   Form1.StringGrid1.Cells[8,i+3]:=FloatToStr(exp(x[i]));
  Form1.StringGrid2.Cells[1,1]:=FloatToStr(exp(x[0])+exp(x[3])*3+exp(x[4])*2+exp(x[5])+exp(x[6])*2);
  Form1.StringGrid2.Cells[1,2]:=FloatToStr(exp(x[0])*2+exp(x[1])+exp(x[2])*2+exp(x[3])+exp(x[5]));
  Form1.StringGrid2.Cells[1,3]:=FloatToStr(exp(x[3])+exp(x[6]));
  Form1.StringGrid2.Cells[1,4]:=FloatToStr(exp(x[1])-exp(x[3])-exp(x[5]));
  Form1.StringGrid2.Cells[1,5]:=FloatToStr(-log10(exp(x[1])));
  Form1.StringGrid2.Cells[1,6]:=FloatToStr(GibbsDuhem(x,G0,n)*R*T);
  Form1.StringGrid2.Cells[1,7]:=FloatToStr(sqrt(2*P_new/rk));
  Form1.StringGrid2.Cells[2,1]:=FloatToStr(g[0]);
  Form1.StringGrid2.Cells[2,2]:=FloatToStr(g[1]);
  Form1.StringGrid2.Cells[2,3]:=FloatToStr(g[2]);
  Form1.StringGrid2.Cells[2,4]:=FloatToStr(g[3]);
  Form1.StringGrid2.Cells[2,5]:=FloatToStr(g[4]);
  Form1.StringGrid2.Cells[2,6]:=FloatToStr(g[5]);
  Form1.StringGrid2.Cells[2,7]:=FloatToStr(g[6]);
  Form1.StringGrid2.Cells[3,1]:=FloatToStr(Lag[0]);
  Form1.StringGrid2.Cells[3,2]:=FloatToStr(Lag[1]);
  Form1.StringGrid2.Cells[3,3]:=FloatToStr(Lag[2]);
  Form1.StringGrid2.Cells[3,4]:=FloatToStr(Lag[3]);
  Form1.StringGrid2.Cells[3,5]:=FloatToStr(Lag[4]);
  TimeEdit.Text:=FormatDateTime('hh:nn:ss:zzz', Now-TimeStart);         //Calculate the duration of calculation
 except
  Application.MessageBox('Fatal error #0001 happened!  ','Lagrange Multipliers',MB_OK);
  Abort;
 end;
end;

procedure TForm1.FormClose(Sender: TObject; var Action: TCloseAction);
begin
 WriteIni;
end;

procedure TForm1.FormPaint(Sender: TObject);
begin
 ReadIni;
end;

end.

[мат-ламер]

alenov
Ваш текст программы я точно не осилю. Во-первых, тщательно протестируйте программу на простых примерах. Допустим, простая квадратичная функция небольшой размерности минимизируется методом сопряжённых градиентов за число шагов, не превышающее размерность пространства. (Считаем, что функция простая и ошибки не накапливаются). Во-вторых. Допустим используется метод модифицированной функции Лагранжа. Там двойная итерация. На внешней итерации пересчитываются множители Лагранжа. На внутренней итерации минимизируется собственно сама эта мод. функция Лагранжа. Так эту внутреннюю минимизацию не обязательно выполнять точно. Что может сильно сократить вычисления.

[alenov]

Спасибо! Попробую простые примеры.

[B@R5uk]

alenov в сообщении #776599 писал(а):

Сконструировал вот такую целевую функцию:

$$ \min P(x^{(k)}, r^{(k)})=\sum(x^{(k)}_i(G_0_i+RT\ln(x^{(k)}_i/\sum(x^{(k)}_i))))+0.5r^{(k)} \sum(h_j(x^{(k)})^2)-r^{(k)} \sum(\ln(-x^{(k)}_i)) $$

Что-то берут меня сомнения, что такую функцию можно минимизировать численно. Лагранж придумал свой метод для поиска условного экстремума в те времена, когда численная математика считалась если не чем-то еретическим, то дурным тоном точно. Его метод создан для аналитического вычисления условного экстремума функции. Он совсем не подходит для численного вычисления экстремума вот в такой простой форме как у вас (покоординатный спуск). И, мне кажется, ни для какого другого численного метода он тоже не подходит, во всяком случае без каких-либо существенных доработок, меняющих его суть.

Почему, я так считаю? Предположим программа находит значения чисел $x^{(k)}, r^{(k)}$, такие, что ваша целевая функция $P$ минимальна. Какие тогда эти числа? Если $x^{(k)}$ доставляют минимум исходной целевой функции $G$ и удовлетворяют наложенным на них дополнительным условиям, то значения выражений $h_j(x^{(k)})$ равны нулю. Следовательно значения $r^{(k)}$ могут быть любыми. Но функция, минимизирующая целевую функцию $P$ возвращает какие-то конкретные $r^{(k)}$ (вы же нигде не прорабатывали такую возможность, что эти множители могут быть любыми?), следовательно, скорее всего функция не правильно вычисляет значения $x^{(k)}$. Один из неправильных решений может быть тот случай, когда один или несколько из множителей Лагранжа равен в результате нулю, а соответствующее ему условие на параметры $x^{(k)}$ не выполняется. Если при этом все условия не выполняются, то $x^{(k)}$ указывают глобальный экстремум функции $G$. Это самый вероятный случай, на мой взгляд, появления неправильного результата.

alenov в сообщении #776599 писал(а):

Значение r у меня пока всегда единица.

Вот важный момент. В этом случае значения $x^{(k)}$ не будут соответствовать ни глобальному экстремуму функции $G$, ни условному.

Ещё мне не понятно, зачем у функции $P$ вот это слагаемое:

alenov в сообщении #776599 писал(а):

$$-r^{(k)} \sum(\ln(-x^{(k)}_i)) $$

Ведь, если величины $x^{(k)}_i$ -- количество вещества, то под логарифмом оказываются отрицательные числа.

[alenov]

B@R5uk в сообщении #812486 писал(а):

alenov в сообщении #776599 писал(а):

Сконструировал вот такую целевую функцию:

$$ \min P(x^{(k)}, r^{(k)})=\sum(x^{(k)}_i(G_0_i+RT\ln(x^{(k)}_i/\sum(x^{(k)}_i))))+0.5r^{(k)} \sum(h_j(x^{(k)})^2)-r^{(k)} \sum(\ln(-x^{(k)}_i)) $$

Это штрафная функция в общем виде. В итоговом варианте я использую не ее. Что я использую, лучше всего посмотреть в коде (хотя это не окончательный вариант, но идеалогически ничего не поменялось). Штрафная функция устроена таким образом, что на начальных иттерациях она маленькая. Алгоритм в основном минимизирует функцию Гиббса. Потом она возрастает и становится намного больше любого значения функции Гиббса. Т.е. сначала находится какое-то приближенное решение, а потом оно подгоняется под ограничивающие условия. Достигается это увеличением числа r. В итоговом варианте у меня сначала r $= 1$, и при удовлетворении определенных условий множится на 2. Есть другие стратегии работы с этим множителем. Часть с минусом у x задействуется только если х отрицательное. В итоговом варианте я это не использую, потому что х всегда положительные.