Приближенное решение СЛАУ с ограничениями (LAPACK || GSL)

[Gcloud]

Здравствуйте!
Сами мы не местные , не математики т.е.
Я программист, но серьезной математикой заниматься особо не приходилось.
Столкнулся я с необходимостью решения системы линейных уравнений; к этому свелась задача. Попробую сформулировать...

Имеется система уравнений
A1.1*X1 + A1.2*X2 +...+ A1.n*Xn = B1
A2.1*X1 + A2.2*X2 +...+ A2.n*Xn = B2
Am.1*X1 + Am.2*X2 +...+ Am.n*Xn = Bm

Для каждого уравнения есть свои граничные условия: p1.1 >= b <= p2.1 ... p1.m >= b <= p2.m . Значения p зависят от B и задаются в виде % от него. Пользователь задает их сам - погрешность(+/-) для каждого уравнения в процентах от B.

И условие: X(1...n) >= 0
Все значения в матрицах A и B, так же >= 0.
Матрицы могут иметь размерности: 25 >= n >= 1, 35 >= m >= 1.
Но наиболее частая ситуация n = 3...5, m = 10...15. Хотя будут и совсем простые, типа n=1...2, m=1...3.
Причем матрицы, в основном, будут иметь вид:
|_26.04___43.821___0______|__103.21
|_1.69____3_______14.654__|__55.12
|_55.80___12.47____54.324_|__101.4
|_0_______31.586___0_____|__92.6
|_4.6_____0.43_____18.221_|__67
|_12.04___0________33.645_|__71
|_0.021___0.72_____0______|__1.2
|_0.006___0.14_____0______|__0.2
|_0.007___0.21_____0______|__0.51
|_0_______0.071____0______|__0.14
|_0_______0.076____0______|__0.03
|_0_______0.004____0______|__0.01

В основном более разреженные снизу. Основная плотность будет в верхней части (m=1...6).

Необходимо найти наиболее близкие решения, т.к. точные решения будут возможны лишь в 1-5% случаев.
Т.е. получаются уравнения вида Am.1*X1 + Am.2*X2 +...+ Am.n*Xn - (bm + [0...p2 || 0...-p1]) => 0 (не знаю как это изобразить математически). Разность левой и правой части с соответствующими допустимыми отклонениями, стремится к нулю.

Разбирался с разными методами решения (Якоби, Гаусса, Йордана-Гаусса, Гаусса-Зейделя, Чолески, QRD, LUD и SVD). Почитав литературу и форум я понял, что для разного вида систем необходимо использовать разные методы.
И вот с вопросом "какие и когда использовать?", разобраться не могу .
Сейчас выбираю каким пакетом воспользоваться - CLAPACK или GSL, еще не определился...

В общем хочу разработать алгоритм решения такого вида уравнений, но моих знаний явно не хватает .

Какие должны быть шаги по определению метода для решения подобных СЛАУ? Как вообще решать системы уравнений с диапазонами вместо точных значений в матрице B?

[Gcloud]

Я так понимаю, нужно внести параметры в правую часть уравнения (B).
Можно ли ее представить в виде системы уравнений Am.1*X1 + Am.2*X2 +...+ Am.n*Xn - ( Bm + p2*X(n+1) - Bm + p1*X(n+2) ) = 0? Каким образом внести эти параметры в систему уравнений, вместе с условием X(1...n) >= 0?

[Yuri Gendelman]

Ваша система - не СЛАУ, а система $2 \cdot m + n$ линейных неравенств с $n$ неизвестными.
Решать ее нужно методами линейного программирования.
Предварительно Вы должны определить критерий оптимизации. Например, можно минимизировать отклонения значений $b_i$ (сумму квадратов или максимум).

Если бы у Вас не было ограничений на $x_i$ а значения $b_i$ были точными, то у Вас была бы система $m$ линейных уравнений с $n$ неизвестными. Т.к. матрица системы - не квадратная, то искать можно было бы только псевдорешение, например, методом SVD.

[Gcloud]

Yuri Gendelman, спасибо. Теперь хоть ясно, что действовать нужно иначе.

Цитата:

Ваша система - не СЛАУ, а система $2 \cdot m + n$ линейных неравенств с $n$ неизвестными.

Честно говоря, не очень понимаю, что Вы имеете ввиду . Почему именно $2 \cdot m + n$? Уравнений может быть как больше, так и меньше кол-ва неизвестных.

Цитата:

Решать ее нужно методами линейного программирования.

А методы решения СЛАУ к ним не относятся ?

Цитата:

Предварительно Вы должны определить критерий оптимизации. Например, можно минимизировать отклонения значений $b_i$ (сумму квадратов или максимум).

Не очень понимаю, как это сделать, ведь минимизировать отклонение конкретного значения $b$ нужно учитывая все неизвестные в $A$ и остальные значения $b$.
Не могли бы Вы дать ссылочку на пример решения подобных неравенств? И возможно ли их решение с использованием CLAPACK или GSL? Если нет, то какую Сишную библиотеку порекомендуете использовать?

Цитата:

Если бы у Вас не было ограничений на $x_i$ а значения $b_i$ были точными, то у Вас была бы система $m$ линейных уравнений с $n$ неизвестными. Т.к. матрица системы - не квадратная, то искать можно было бы только псевдорешение, например, методом SVD.

Насчет $x_i \geqslant 0$, мне кажется я где-то видел алгоритм решения, на основе то ли Гаусса, то ли Якоби. Не помню, много инфы пришлось перерыть по теме.
Я с помощью SVD и пытался решать эти уравнения, но без учета параметров. Получалась полная ерунда на некоторых вариантах. Уравнения очень похожи(размеры матрицы одинаковые и распределение коэффициентов приблизительно такое же, но одно решается более-менее нормально, а второе вообще не слишком логично. И во многом проблема в отсутствии ограничений при расчете. Вот и думал, как бы их внести в уравнение. Видимо никак...

[Yuri Gendelman]

Gcloud в сообщении #460820 писал(а):

Честно говоря, не очень понимаю, что Вы имеете ввиду . Почему именно $2 \cdot m + n$? Уравнений может быть как больше, так и меньше кол-ва неизвестных.

У Вас $m$ "уравнений" вида
$A_{j,1} X_1 + A_{j,2} X_2 + ... + A_{j,N} X_N = B_j$.

Но, поскольку $ Bmin_j \leq B_j \leq Bmax_j$, каждое "уравнение" - это на самом деле - пара неравенств
$A_{j,1} X_1 + A_{j,2} X_2 + ... + A_{j,N} X_N \leq Bmax_j$.
$A_{j,1} X_1 + A_{j,2} X_2 + ... + A_{j,N} X_N \geq Bmin_j$.

Кроме того, у Вас еще есть $n$ условий $ X_i \geq 0$

Gcloud в сообщении #460820 писал(а):

Цитата:

Решать ее нужно методами линейного программирования.

А методы решения СЛАУ к ним не относятся ?

Наиболее известный метод решения задач линейного программирования - симплексный метод. Он широко применяется в задачах экономического планирования. На эту тему есть много литературы.
О конкретных пакетах ничего сказать не могу, я за ними уже давно не слежу.

[Gcloud]

Спасибо за доступное объяснение! Теперь ясно... в среднем система из 27-28 линейных неравенств вида:
$A_{j,1} X_1 + A_{j,2} X_2 + ... + A_{j,N} X_N \leq Bmax_j$.
$A_{j,1} X_1 + A_{j,2} X_2 + ... + A_{j,N} X_N \geq Bmin_j$.
...
$ X_1 \geq 0$
...
$ X_i \geq 0$

Уже ищу инфу по симплекс-методу и вообще методам решения неравенств. Но литературы действительно много.
Если бы кто накидал в тему ссылочек по практической стороне решения систем линейных неравенств, а еще бы ткнул носом в какой библиотеке (C, C++, Java, JS, PHP, Perl) искать средства решения... а то опять зароюсь в теорию в поисках методов решения .

[Gcloud]

Цитата:

Предварительно Вы должны определить критерий оптимизации. Например, можно минимизировать отклонения значений $b_i$ (сумму квадратов или максимум).

Yuri Gendelman, я так понимаю, нужно определить целевую ф-цию... Но я в упор не понимаю как это сделать. У меня задача очень похожа на классическую задачу о диете, только без стоимости продуктов и ее минимизации. Т.е. имеется набор продуктов с содержанием питательных элементов, и есть таблица с нормами пит. элементов. Задача - подобрать оптимальные соотношения масс выбранного набора продуктов по выбранной программе питания.

______________--------------- $A$ ------------________---$B$--_--$p_1, p_2$---___
______________ Хлеб____Молоко___Сыр________Программа питания (погрешность)_
Пит. ценность____120______180_____250____________220 (-10/+9 %)
Жиры___________5________21______51_____________96 (-7/+4 %)
Белки___________7________12______42_____________62 (-5/+3 %)
Углеводы________15_______17______14_____________43 (-5/+5 %)
Витамин С________0,1______0,2______0,32___________0,85 (-15/+12 %)
Витамин Е________0,01_____0,03_____0,04___________0,2 (-35/+45 %)
___________________________________________________________
Вес продукта:_____$x_1$ г____$x_2$ г____$x_3$ г

$p_{1,i} = b_i - p_1 \%$
$p_{2,i} = b_i + p_2 \%$

Подскажите пожалуйста, какую тут можно задать целевую ф-цию?

[Yuri Gendelman]

Gcloud в сообщении #461736 писал(а):

Задача - подобрать оптимальные соотношения масс выбранного набора продуктов по выбранной программе питания.

Если Ваша задача - "подобрать оптимальные соотношения масс", то у Вас уже должен быть критерий оптимальности: какие соотношения масс "более оптимальны, чем другие".
Например, Вы можете минимизировать стоимость продуктов, отклонение от идеальной программы питания, и т.д.

[Gcloud]

Yuri Gendelman, я как раз нашел решение, но все равно спасибо.
И даже нашел замечательный бесплатный мат. пакет линейного программирования - LP_Solve. Как раз он очень удобен для моих целей, он имеет интерфейс для PHP, причем в виде готового PHP-модуля - врапера для библиотеки. Не придется писать самому :).
Еще одно из преимуществ - имеет свою IDE, простую, но удобную, что упрощает обучение пакету и процессы оптимизации расчетов. Не нужно юзать консольные утилиты или компилить код.
Осталось переписать математическое представление задачи в функциональное.

В итоге получилось так:

$x_1 \geqslant 0$
...
$x_n \geqslant 0$

$u_{1 (minimal)} \leqslant b_1 \leqslant u_{1 (maximal)}$
...
$u_{m (minimal)} \leqslant b_m \leqslant u_{m (maximal)}$

$a_{1,1} x_1 + a_{1,2} x_2 + ... + a_{1,n} x_n = b_1$
...
$a_{m,1} x_1 + a_{m,2} x_2 + ... + a_{m,n} x_n = b_m$

|c_1-b_1 + c_2-b_2 + c_3-b_3 + ... + c_n-b_n| \rightarrow \min$

Где $c_1...c_n$ - идеальные значения, взятые из программы питания
а $u_{m (minimal)}$, $u_{m (maximal)}$ - граничные значения для $c_n$, заданные в программе питания.
Уравнения сортируются по значениям $b$ в порядке возрастания.

Решает очень хорошо, и, что немаловажно, быстро.

Т.к. наличие решения необходимо "любой ценой", то при неудачной попытке решения, нужно как-то находить уравнение, дающее наибольшее отклонение $c_n - b_n$, исключать его из условия, и снова пытаться найти решение. И так до тех пор, пока не будет найдено решение.
Вот только непонятно, можно ли (и как?) по данным из недорешенного уравнения, определить "проблемное" уравнение?

[Yuri Gendelman]

Ваше сообщение о решателе задач линейного программирования я скопировал в раздел "Интернет-ресурсы" (Тема "Численные методы").

Gcloud в сообщении #462014 писал(а):

Вот только непонятно, можно ли (и как?) по данным из недорешенного уравнения, определить "проблемное" уравнение?

Не знаю, возможно ли это непосредственно. Это зависит от структуры пакета.

Но Вы можете:
- увеличить диапазоны для всех уравнений и решить новую задачу,
- найти уравнение с максимальным отклонением и исключить его,
- восстановить старые диапазоны и решить новую задачу.

Этот цикл можно повторить несколько раз.

[Gcloud]

Цитата:

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

Дело в том, что расширив диапазон даже одного уравнения, пересчет может пойти совсем по другим путям, и вместо использования всех продуктов в списке, некоторые могут быть исключены из него, а некоторые устремятся за пределы, которые реально не должны быть доступны, что приведет к проблемам на шаге №3. Выяснил при тестировании модели. Я все таки надеюсь на то, что библиотека позволяет произвести какую-то оценку по результатам расчетов или каким то отдельным методом. В результатах есть такие данные: "Objective", "Constraints", "Sensitivity". Но при неудачной попытке решения, они пустые. Быть может есть какой-то математический прием для оценки? Если таковой имеется, то возможно он реализован в библиотеке.

[Yuri Gendelman]

Gcloud в сообщении #462207 писал(а):

пересчет может пойти совсем по другим путям,...

Да, задача линейного программирования может оказаться некорректной (нет непрерывной зависимости решения от исходных данных на всей области параметров). Пробуйте малые изменения диапазонов.

Чтобы продукты не исключались, замените условия $X_j \geq 0$ на $X_j > 0$.

Gcloud в сообщении #462207 писал(а):

Быть может есть какой-то математический прием для оценки?

Если в момент обнаружения ошибки пакет сохраняет последние вычисленные значения $b_i, X_j$, попытайтесь их использовать.

[Gcloud]

Цитата:

Пробуйте малые изменения диапазонов.
Чтобы продукты не исключались, замените условия $X_j \geq 0$ на $X_j > 0$.

Спасибо, буду пробовать. Написал скриптик для анализа работы алгоритма и для поиска проблемных уравнений, использовал метод релаксации:
$a_{1,1} x_1 + a_{1,2} x_2 + ... + a_{1,n} x_n + e_1= b_1$
$a_{m,1} x_1 + a_{m,2} x_2 + ... + a_{m,n} x_n + e_m = b_m$

$|c_1-b_1 + c_2-b_2 + c_3-b_3 + ... + c_n-b_n + k e_1 + ... k e_m| \rightarrow \min$
Где $k$ - достаточно большое число, больше любого члена в уравнениях и неравенствах у меня $k = 1000$.
В результате получаю значения $e$ для каждого уравнения, которые являются индикаторами отклонений, необходимых для коррекции верхних или нижних границ. Если значение отлично от 0, то соответствующее уравнение нуждается в корректировке границ. Правда алгоритм предлагает не самые оптимальные пути решения - находит результаты с достаточно большими отклонениями, разрешив которые решение будет весьма неудовлетворительным. Более оптимальные варианты решения существуют, проверял подбором вручную. Вот как бы этот подбор алгоритмизировать, да чтоб не перебрать с итерациями...
Есть у меня подозрения, что алгоритм решения в самой библиотеке выбирается не оптимально, я в экселе за 5 минут подбираю гораздо более оптимальные решения. Быть может нейросеть обучить?

Т.к. я не силен в теории, то не знаю, какие можно применять еще методы. А ведь они должны быть, ведь намного более сложные задачи оптимизации решаются во множестве софта, с использованием все тех же библиотек (lp_solve).
Вот только документация у lp_solve не такая подробная, как хотелось бы. Например так и не нашел как через API ввести в уравнение неизвестную, являющуюся неизвестным в неравенстве. Т.е. значения $b$ в моем примере. Пришлось написать генератор мат.модели, который строит ее на основе двумерного массива и сохраняет ее в файл, а затем библиотека его открывает и решает.

Через debug-режим можно писать в файл такую инфу:

Цитата:

duals dualsfrom dualstill R
0 -0.002 0.0016666666666667 -0.0036666666666667
0 -1.0E+30 1.0E+30 -2.0E+30
0 -1.0E+30 1.0E+30 -2.0E+30
0 -0.075 1.0E+30 -1.0E+30
0 -1.0E+30 1.0E+30 -2.0E+30
0 -1.0E+30 1.0E+30 -2.0E+30
0 -34.6 6.7177254472435 -41.317725447244

Но как ее использовать я не понимаю .

[Yuri Gendelman]

Gcloud в сообщении #463161 писал(а):

$|c_1-b_1 + c_2-b_2 + c_3-b_3 + ... + c_n-b_n + k e_1 + ... k e_m| \rightarrow \min$

Я бы выбрал такой критерий (минимизировать максимальное отклонение):
$\max_{1\leq j \leq m}{(\max{(|b_j-bmax_j|, |b_j-bmin_j|)})} \rightarrow \min$

Gcloud в сообщении #463161 писал(а):

Например так и не нашел как через API ввести в уравнение неизвестную, являющуюся неизвестным в неравенстве. Т.е. значения $b$ в моем примере.

Эта фраза мне не понятна.
Обычно на входе задается начальное приближение ($x_{i}^{0}, b_{j}^{0}$), а на выходе получается решение ($x_i, b_j$).