Текст подпрограммы и версий
mnp1r_p.zip 		Тексты тестовых примеров
tmnp1r_p.zip

Подпрограмма:  MNP1R (модуль MNP1R_p)

Назначение

Решение задачи квадратичного программирования при наличии линейных ограничений на переменные методом покоординатного спуска.

Математическое описание

Решается задача квадратичного программирования 

     min  { (p, x) + (Cx, x)  |  Ax ≤ b } ,

где  p и  x - векторы длины  n,  C - симметричная положительно определенная матрица размерности  n * n,  A - матрица размерности  m * n,  b - вектоp длины  m.

Матрица  C задается в компактной форме, т.е. представляется в виде вектоpа длины  n (n + 1)/2, который состоит из элементов нижнего треугольника матрицы, выписанных последовательно по строкам.

Для получения решения исходной задачи строится двойственная ей задача: найти 

     min  { φ(u) = (h, u) + (Gu, u)  |  u ≥ 0 } ,

где  h = (1/2)AC - 1 p + b,  G = (1/4)AC - 1 A*.

Для решения двойственной задачи используется метод покоординатного спуска. Для одномерной минимизации функции  φ (u) вдоль направления спуска используется метод квадратичной аппроксимации.

Некоторая вычисленная точка  uk ≥ 0 считается точкой минимума функции  φ (u), если выполнено хотя бы одно из следующих условий:

1.  | uik - uik - 1 | ≤ XE для всех  i =1, ..., m,  где  uk = (u1k, ..., umk) - точка, полученная на  k - ой итерации метода, а XE - заданная точность решения задачи по аргументу;
2.  | φ (uk) - φ (uk - 1) | ≤ FE, где  uk - точка, вычисленная на  k - ой итерации метода, а FE - заданная точность решения задачи по функционалу.

Решение исходной задачи вычисляется по формуле 

     x(u) = - (1/2) C -1 (A*u + p) .

Кюнци Г.П., Крелле B., Нелинейное программирование, Изд - во "Советское радио", M., 1965.

Использование 

procedure MNP1R(var M :Integer; var N :Integer; var A :Array of Real;
                var B :Array of Real; var U :Array of Real;
                var P :Array of Real; FSTEP :Real; IPAR :Integer;
                var MAXK :Integer; XE :Real; var FE :Real;
                var KOUNT :Integer; var CO :Array of Real;
                var G :Array of Real; var I0 :Array of Integer;
                var IERR :Integer);


   Параметры
M - число стpок матрицы  A (тип: целый);
N - число столбцов матрицы  A (тип: целый);
A - вещественный двумерный массив размерности  M * N;
B - вещественный вектоp длины  M;
U - вещественный вектоp длины  max {6M + N + 11; N (N + 1)/2 }; на входе первые N (N + 1)/2 компонент вектоpа содержат компактную запись матрицы  C;
P - вещественный вектоp длины  N; на входе содержит компоненты вектоpа  p, а на выходе содержит компоненты решения исходной задачи;
FSTEZ - начальная длина шага (тип: вещественный);
IPAR - параметр, управляющий вариантом покоординатного спуска: если IPAR = 1, то используется вариант с ускоренным дроблением шага (тип: целый);
MAXK - целая переменная, на входе равная максимально допустимому числу вычислений функции, а на выходе - числу фактически выполненных вычислений функций;
XE - заданная точность решения задачи по аргументу (тип: вещественный);
FE - заданная точность решения задачи по функционалу (тип: вещественный);
KOUNT - целая переменная, содержащая число фактически выполненных итераций метода;
CO - вещественный вектоp длины N (N + 1)/2, используемый как рабочий;
G - вещественный вектоp длины M (M + 1)/2, используемый как рабочий;
I0 - целочисленный вектоp длины  M, используемый как рабочий;
IERR - целая переменная, указывающая причину окончания процесса:
IERR= 0 - если найден минимум функции с заданной точностью по аргументу или по функционалу;
IERR= 4 - если выполнено заданное максимальное число вычислений функции, но точность не была достигнута;
IERR=65 - если матрица  C не является положительно определенной. 
   Версии:  нет 
   Вызываемые подпрограммы
AIH1R - подпрограмма обращения положительно определенной симметричной матрицы с компактной формой представления методом квадратного корня. 
   Замечания по использованию
 

Матрица  C должна быть положительно определенной.

Используются служебные подпрограммы: MNP11, MNP12, MNP13, MNP14, MNP15, MNP16.

Пример использования 

   Найти
         min  [ Q(x) = 0.5x12 + 0.5x22 - x1 - 2x2 ]

   при ограничениях

         2x1 + 3x2 ≤ 6 , 
           x1 + 4x2 ≤ 5 , 
           x1          ≥ 0 , 
                    x2 ≥ 0 .

   Точка минимума    x*  =  (13/17, 18/17) .


Unit TMNP1R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, MNP1R_p;

function TMNP1R: String;

implementation

function TMNP1R: String;
var
KOUNT,IERR,_i :Integer;
IO :Array [0..3] of Integer;
СО :Array [0..2] of Real;
G :Array [0..9] of Real;
const
ХЕ :Real = 1.0E-8;
FE :Real = 1.0E-9;
M :Integer = 4;
N :Integer = 2;
FSTEZ :Real = 1.0;
IPAR :Integer = 1;
МАХК :Integer = 250;
A :Array [0..7] of Real = ( 2.0,1.0,-1.0,0.0,3.0,4.0,0.0,-1.0 );
B :Array [0..3] of Real = ( 6.0,5.0,0.0,0.0 );
P :Array [0..1] of Real = ( -1.0,-2.0 );
U :Array [0..36] of Real = ( 0.5,0.0,0.5,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0 );
begin
Result := '';  { результат функции }
MNP1R(M,N,A,B,U,P,FSTEZ,IPAR,MAXK,
     XE,FE,KOUNT,CO,G,IO,IERR);
Result := Result + Format('%s',['  IERR=']);
Result := Result + Format('%2d ',[IERR]);
Result := Result + Format('%s',[' KOUNT=']);
Result := Result + Format('%3d ',[KOUNT]);
Result := Result + Format('%s',['  MAXK=']);
Result := Result + Format('%4d ',[MAXK]);
Result := Result + Format('%s',['      P(1)=']);
Result := Result + #$0D#$0A;
for _i:=0 to 1 do
 begin
  Result := Result + Format('%20.16f ',[P[_i]]);
  if ( ((_i+1) mod 2)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TMNP1R',Result);  { вывод результатов в файл TMNP1R.res }
exit;
end;

end.

Результаты:

      IERR  =  0

      KOUNT  =  95
      MAXK  =  117

      P(1)  =  7.646729 - 01
      P(2)  =  1.058691 + 00