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

Подпрограмма:  mnp1r_c

Назначение

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

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

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

     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.

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

    int mnp1r_c (integer *m, integer *n, real *a, real *b, real *u,
             real *p, real *fstep, integer *ipar, integer *maxk, real *xe,
             real *fe, integer *kount, real *co, real *g, integer *i0,
             integer *ierr)

   Параметры
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, а на выходе содержит компоненты решения исходной задачи;
fstep - начальная длина шага (тип: вещественный);
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 - подпрограмма обращения положительно определенной симметричной матрицы с компактной формой представления методом квадратного корня. 
   Замечания по использованию
 

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

Используются служебные подпрограммы: mnp11_c, mnp12_c, mnp13_c, mnp14_c, mnp15_c, mnp16_c.

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

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

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

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

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

int main(void)
{
    /* Initialized data */
    static float xe = 1e-8f;
    static float p[2] = { -1.f,-2.f };
    static float u[37] = { .5f,0.f,.5f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
          0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
          0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f };
    static float fe = 1e-9f;
    static int m = 4;
    static int n = 2;
    static float fstep = 1.f;
    static int ipar = 1;
    static int maxk = 250;
    static float a[8] /* was [4][2] */ = { 2.f,1.f,-1.f,0.f,3.f,4.f,0.f,-1.f };
    static float b[4] = { 6.f,5.f,0.f,0.f };

    /* Local variables */
    static int ierr;
    extern int mnp1r_c(int *, int *, float *, float *, float *, float *,
                       float *, int *, int *, float *, float *, int *,
                       float *, float *, int *, int *);
    static float g[10];
    static int kount;
    static float co[3];
    static int io[4];

    mnp1r_c(&m, &n, a, b, u, p, &fstep, &ipar, &maxk, &xe, &fe, &kount, co,
            g, io, &ierr);

    printf("\n %5i \n", ierr);
    printf("\n %5i \n", kount);
    printf("\n %5i \n", maxk);
    printf("\n %16.7e %16.7e \n", p[0], p[1]);
    return 0;
} /* main */


Результаты:

      ierr  =  0

      kount  =  95
      maxk  =  117

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