алгоритм оптимизации

Заранее известны матрицы A, B, CM, DM, вектора x, s, s*, cm, dm, u0, все элементы неотрицательные, число epsilon для точности. u-вектор с неизвестными элементами известной размерности p. 

CM-матрица правых частей, cm-вектор левых частей ограничения CM*(A*x+B*u-s)<=cv (1) . 

DM-матрица правых частей, dm-вектор левых частей ограничения DM*u<=dv (2) . 

Из исходных данных строится функция, зависящая от u, евклидовой нормы вектора A*x+B*u-s-s*, а также градиент этой функции. 

1)Для точки u0 выбираются активные и неактивные ограничения из (1) и (2) . 

2)Находится вектор направления d решением задачи линейного математического программирования симплекс-методом, где целевая функция - найденный ранее градиент в точке u0 умноженный на вектор d(d имеет размерность градиента), а ограничения строятся, как матрица правых частей активных ограничений * d<= вектор левых частей активных ограничений, а также дополнительное ограничение 0<=d<=1

3)Если в найденном d целевая функция по модулю меньше epsilon, алгоритм заканчивается, выводится u0, xt=A*x+B*u0-s-s*, значение функции, заданной евклидовой нормой выше в точке u0. 

4)В противном случае выбирается новая точка uk по правилу uk=u0+l*d, где l - число. Строится функция F(l), полученная подстановкой uk в функцию, найденной ранее через евклидову норму. l ищется методом равномерного поиска минимизации F(l) с ограничением 0<=l<=l_max. l_max=min(bi/mi), где bi-элементы вектора b: b=вектор правых частей неактивных ограничений-матрица левых частей неактивных ограничений*u0, mi-элементы вектора m: m=матрица левых частей неактивных ограничений*d.

5) вместо точки u0 берëм uk и повторяем алгоритм с шага 1 (переопределяем активные и неактивные ограничения, вектор направления, при необходимости ищем новую точку) пока в каком-то d целевая функция по модулю не будет меньше epsilon

нужно реализовать решение данного алгоритма на языке python