Автор: Евгений Шлыков
Группа: 596
Семинарист: Алексей Глибичук
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
EPS = 1e-6
MAXITER = 100000
Реализация симплекс-метода немного отличается от того, как нас учили, поскольку я опирался на книгу P. R. Thie, G. H. Keough "An Introduction to Linear Programming and Game Theory".
Написать симплекс-метод и поддержку следующей функциональности:
n = 2;draw=True и n = 2 выводится график из множества точек. Синие точки обходились первыми, красные в середине пути, желтые последними. Иными словами, порядок - как в тепловизоре при повышении температуре.n и m симплекс-метод работает за O(1), что и требовалось по условию.Функция, которая реализовывает метод Бланда и лексикографический метод для выбора покидающей и входящей переменной.
Здесь стоит отметить, что если входящую переменную можно найти, а выходящую нет, то значит, задача не ограничена. А если нельзя найти даже входящую, то тогда мы нашли оптимальное решение.
In [2]:
def find_entering_leaving(st, phase, method='blend'):
"""
Находит пару входящей и покидающей переменных.
"""
k = 3 - phase
# Отбираем переменные, которые можно ввести в базис.
mo = np.ma.masked_where(st[-1, :-1] <= EPS, st[-1, :-1], copy=False)
if mo.count() == 0:
# Задача недопустимая
return None, None
# Выбираем входящую переменную
if method == 'blend':
ent = np.where(np.logical_not(mo.mask))[0][0]
else:
ent = np.argmax(mo)
ma = np.ma.masked_where(st[:-k, ent] <= EPS, st[:-k, ent], copy=False)
if ma.count() == 0:
# Задача неограниченная
return ent, None
mb = np.ma.masked_where(st[:-k, ent] <= EPS, st[:-k, -1], copy=False)
mc = mb / ma
# Выбираем покидающую переменную
if method == 'blend':
lea = np.argmin(mc)
else:
rows = np.empty((0, st.shape[1] + 1), float)
mm = np.ma.masked_where(mc != np.min(mc), mc, copy=False)
for lea in np.ma.where(np.logical_not(mm.mask))[0]:
row = np.array([st[lea, :] / st[lea][ent]])
row = np.append(row, [[lea]], axis=1)
rows = np.append(rows, row, axis=0)
# Сортируем строчки матрицы в лексикографичесом порядке
# и возвращаем номер наименьшей. Затем получаем номера переменных.
ind = np.lexsort(rows[:, ::-1].T)[0]
lea = int(rows[ind][-1])
return ent, lea
Функция, которая по данной симплекс-таблице высчитывает координаты соответствующей точки.
In [3]:
def add_point(points, st, basis, nvars, k):
"""
По имеющейся симплекс-таблице вычислить текущую базисную точку.
"""
point = np.zeros(nvars, dtype=np.float64)
point[basis[:k]] = st[:k, -1]
return np.append(points, np.array([point]), axis=0)
Функция, которая получает симплекс-таблицу и совершает все итерации для вычисления оптимального решения.
In [4]:
def simplex_method(st, basis, nvars, phase=2, method='blend'):
"""
Применяет симплекс-метод к данной задаче линейного программирования
в стандартной форме.
"""
iters = 0
points = np.zeros(shape=(0, nvars))
m, n = st.shape
k = m - 3 + phase
if phase == 1:
points = add_point(points, st, basis, nvars, k)
if phase == 2:
# Если какая-то aritificial-переменная осталась в базисе после фазы 1,
# то ее необходимо оттуда вывести. См. пример Thie 3.7.2.
# Если закомментировать этот кусок кода, то алгоритм не сработает,
# так как в конце останется aritificial-переменная в базисе, которую
# будет невозможно оттуда вывести.
for leave in [row for row in range(basis.size) if basis[row] > n - 2]:
for_enter = [col for col in range(n - 1) if st[leave, col] != 0]
enter = for_enter[0]
st[leave, :] = st[leave, :] / st[leave][enter]
for i in range(m):
if i != leave:
st[i, :] = st[i, :] - st[leave, :] * st[i, enter]
basis[leave] = enter
iters += 1
points = add_point(points, st, basis, nvars, k)
while iters < MAXITER:
enter, leave = find_entering_leaving(st, phase, method)
if (enter is None) or (leave is None):
if enter is not None:
# Задача не ограничена
return None, None
break
st[leave, :] = st[leave, :] / st[leave][enter]
for i in range(m):
if i != leave:
st[i, :] = st[i, :] - st[leave, :] * st[i, enter]
basis[leave] = enter
iters += 1
points = add_point(points, st, basis, nvars, k)
return iters, points
Заглавная функция, именно она вызывается для решения задачи. Она подготовливает сипдекс-таблицы для обеих фаз и использует метод выше, а также обрабатывает возвращенный результат.
In [5]:
def solve_lin_prog(A, b, c, draw=False, method='blend', start_point=None):
"""
Создает симплекс-таблицу и применяет к ней симплекс-метод.
"""
A = np.array(A)
b = np.array(b)
c = np.array(c)
start_point = np.array(start_point) if start_point is not None else None
b = np.ravel(b) # Приводим к нужному типу.
n = len(c) # Количество переменных.
m = len(b) # Количество slack-переменных.
draw = draw and n == 2
if start_point is not None:
# Сделаем замену координат: y[i] = x[i] - start_point[i]
# А затем y[i] = p[i] - q[i], где p[i], q[i] >= 0
b = b - A @ start_point
b = np.append(b, start_point)
c = np.append(c, -c)
E = np.eye(n)
A = np.bmat([[A, -A], [-E, E]])
m += n
n *= 2
k = np.sum(b < 0) # Количество artificial-переменных.
v = n + m + k # Общее число переменных
# Заполняем симплекс-таблицу.
st = np.zeros([m + 2, v + 1]) # Симплекс-таблица.
constr = st[:-2, -1] # Ссылка на столбец ограничений.
st[-2, :n] = c # Целевая функция.
st[-2, -1] = 0 # Текущий ответ.
st[:m, :n] = A # Сама матрица.
constr[:m] = b # Ограничения в таблице.
np.fill_diagonal(st[:m, n:n + m], 1) # Единичная подматрица для базиса.
# Заполняем базис. Если строка содержит отрицательное ограничеие,
# то к ней добавляется artificial-переменная. А сама строка
# умножается на -1, чтобы сделать ограничение положительным.
basis = np.zeros(m, dtype=int) # Базис.
artificialed = np.zeros(k, dtype=int) # Строки с artificial-переменными.
im = 0 # Индекс slack-переменной.
ik = 0 # Индекс artificial-переменной.
for i in range(m):
if constr[i] >= 0:
# В базис входит slack-переменная.
basis[i] = n + im
im += 1
continue
# Artificial-переменная должна попасть в базис.
basis[i] = n + m + ik
artificialed[ik] = i
ik += 1
constr[i] *= -1
st[i, :-1] *= -1
st[i, basis[i]] = 1
st[-1, basis[i]] = -1 # Она входит со знаком минус.
# Суммируем по всем добавочным переменным в последнюю строку -
# целевую функцию двойственной задачи.
for r in artificialed:
st[-1, :] = st[-1, :] + st[r, :]
# Начинаем первую фазу.
nit1, pts1 = simplex_method(st, basis, v, phase=1, method=method)
# Целевая функция должна быть равна 0 после окончания первой фазы.
# Иначе не удалось найти подходящую стартовую базисную точку.
if (abs(st[-1, -1]) >= EPS) or (nit1 is None) or (pts1 is None):
return None, None, None
# Удаляем строку целевой функции для двойственной задачи
# и столбцы aritifical-переменных.
st = st[:-1, :]
st = np.delete(st, np.s_[n + m:n + m + k], 1)
# Начинаем вторую фазы.
nit2, pts2 = simplex_method(st, basis, v, phase=2, method=method)
if (nit2 is None) or (pts2 is None):
return None, None, None
# Нарисуем путь перебора точек.
if draw:
points = np.vstack((pts1, pts2))[:, :n]
if start_point is not None:
points = points[:, :n // 2] - points[:, n //2:] + start_point
color = np.arange(points.shape[0]) / points.shape[0]
plt.figure(figsize=(6, 6))
plt.plot(points[:, 0], points[:, 1], ls='-', c='black', lw=2, zorder=-1)
plt.scatter(points[:, 0], points[:, 1], c=color, s=500, cmap='plasma')
plt.show()
# Выписываем теперь ответ.'
solution = np.zeros(v)
solution[basis[:m]] = st[:m, -1]
x = solution[:n]
if start_point is not None:
x = x[:n // 2] - x[n // 2:] + start_point
c = c[:n // 2]
opt = x @ c
return x, opt, nit1 + nit2
Функция, которая проверяет, что приведенное решение оптимально.
In [6]:
def is_optimal(A, b, c, x):
"""
За константное при фиксированных n, m проверяет оптимальность решения.
"""
c = np.array(c)
x = np.array(x)
return abs(solve_lin_prog(A, b, c)[1] - np.dot(c, x)) <= EPS
Тесты, взятые из книжек и из условия этого практикума. При n = 2 также выводятся пройденные базисные точки. Если ничего не выведено, тест пройден (assert срабатывает, когда условие ложно).
In [7]:
# Тест из условия
_c = [5, 1]
_A = [[1, 2], [2, 0.5]]
_b = [5, 8]
res = solve_lin_prog(_A, _b, _c, draw=True)
print(res)
assert(abs(res[1] - 20) <= EPS)
assert(is_optimal(_A, _b, _c, res[0]))
In [8]:
# Вандербей стр. 17
_c = [-2, -1]
_A = [[-1, 1], [-1, -2], [0, 1]]
_b = [-1, -2, 1]
res = solve_lin_prog(_A, _b, _c, draw=True)
print(res)
assert(abs(res[1] - (-3)) <= EPS)
assert(is_optimal(_A, _b, _c, res[0]))
In [9]:
# Вандербей 2.1
_c = [6, 8, 5, 9]
_A = [[2, 1, 1, 3], [1, 3, 1, 2]]
_b = [5, 3]
res = solve_lin_prog(_A, _b, _c)
print(res)
assert(abs(res[1] -17) <= EPS)
assert(is_optimal(_A, _b, _c, res[0]))
In [10]:
# Вандербей 2.2
_c = [2, 1]
_A = [[2, 1], [2, 3], [4, 1], [1, 5]]
_b = [4, 3, 5, 1]
res = solve_lin_prog(_A, _b, _c, draw=True)
print(res)
assert(abs(res[1] - 2) <= EPS)
assert(is_optimal(_A, _b, _c, res[0]))
In [11]:
# Вандербей 2.3
_c = [2, -6, 0]
_A = [[-1, -1, -1], [2, -1, 1]]
_b = [-2, 1]
res = solve_lin_prog(_A, _b, _c)
print(res)
assert(abs(res[1] - (-3)) <= EPS)
assert(is_optimal(_A, _b, _c, res[0]))
In [12]:
# Вандербей 2.4
_c = [-1, -3, -1]
_A = [[2, -5, 1], [2, -1, 2]]
_b = [-5, 4]
res = solve_lin_prog(_A, _b, _c, draw=True)
print(res)
assert(abs(res[1] - (-3)) <= EPS)
assert(is_optimal(_A, _b, _c, res[0]))
In [13]:
# Вандербей 2.5
_c = [1, 3]
_A = [[-1, -1], [-1, 1], [1, 2]]
_b = [-3, -1, 4]
res = solve_lin_prog(_A, _b, _c, draw=True)
print(res)
assert(abs(res[1] - 5) <= EPS)
In [14]:
# Вандербей 2.8
_c = [3, 2]
_A = [[1, -2], [1, -1], [2, -1], [1, 0], [2, 1], [1, 1], [1, 2], [0, 1]]
_b = [1, 2, 6, 5, 16, 12, 21, 10]
res = solve_lin_prog(_A, _b, _c, draw=True, start_point=[4, 2])
print(res)
assert(abs(res[1] - 28) <= EPS)
assert(is_optimal(_A, _b, _c, res[0]))
In [15]:
# Thie 3.7.2
_c = [-1, -4, -3, -2]
_A = [[1, 2, 0, 1], [2, 1, 1, 0], [-1, 4, -2, 3],
[-1, -2, 0, -1], [-2, -1, -1, 0], [1, -4, 2, -3]]
_b = [20, 10, 40, -20, -10, -40]
res = solve_lin_prog(_A, _b, _c)
assert(abs(res[1] - (-35)) <= EPS)
assert(is_optimal(_A, _b, _c, res[0]))
y = a0 + a1 * t + a2 * sin(t).y_real.y_real надо добавить значения, распределенные нормально N(0, Sigma), получив тем самым y_corr.Sigma.y_corr - y_real.a0, a1, a2, выбрав один из трех методов:L1-норме;L2-норме;LINF-норме.Sigma.y_corr: добавить к первому 50, а из второго отнять 50 - и посмотреть, как это отразится на результате.Функция, которая вычисляет y_real по данным коэффициентам a в точках t.
In [16]:
def get_y_real(a, t):
"""
Вычисляет истинное значение функции.
"""
return a[2] * np.sin(t) + a[1] * t + a[0]
Функция, которая к y_real добавляет искажения, распределенные как N(0, Sigma), возвращая y_corr.
In [17]:
def get_y_corr(y_real, sigma=1.0):
"""
Вычисляет искаженные данные по настоящим.
"""
return y_real + np.random.normal(0, sigma, len(y_real))
Функция, которая рисует графики самой функции, искаженных точек и графики функций, где вместо параметров подставлены посчитанные параметры.
In [18]:
def draw(a_real, a_corr, y_corr, t, maxt, sigma, pm50):
"""
Рисует график исходной кривой и ее приближений тремя методами.
"""
x = np.linspace(0, maxt, 100)
y = get_y_real(a_real, x)
plt.figure(figsize=(18, 6))
plt.title('Sigma = {}, ЕстьВыбросы = {}'.format(sigma, pm50))
plt.scatter(t, y_corr, alpha=0.25, s=25, c='y')
plt.plot(x, y, c='y', label='Функция', lw=2)
plt.plot(x, get_y_real(a_corr[0], x), c='r', label='По L2-норме',)
plt.plot(x, get_y_real(a_corr[1], x), c='g', label='По L1-норме')
plt.plot(x, get_y_real(a_corr[2], x), c='b', label='По LINF-норме')
plt.ylim((-np.min(y) - sigma, np.max(y) + sigma))
plt.legend()
plt.show()
Функция, которая получает оптимальные значения параметров для выбранного метода.
In [19]:
def get_params(y_corr, t, method=0):
"""
По набору точек (t[i], y_corr[i]) находит значения параметров a, используя:
method = 0: корень суммы квадратов значений вектора невязок (L2-норма);
methid = 1: сумма модулей значений вектора невязок (L1-норма);
method = 2: максимальный модуль значения вектора невязок (LINF-норма).
"""
m = len(t)
B = np.matrix([[1, ti, np.sin(ti)] for ti in t])
y_corr = np.matrix(y_corr).T
b = np.bmat([[y_corr], [-y_corr]])
if method == 0:
return ((B.T @ B).I @ B.T @ y_corr).A1
elif method == 1:
E = np.eye(m, dtype=float)
A = np.bmat([[-E, B], [-E, -B]])
c = np.bmat([np.ones(m), [0, 0, 0]]).A1
else:
e = np.ones(len(t), dtype=float).reshape(len(t), 1)
A = np.bmat([[-e, B], [-e, -B]])
c = np.array([1, 0, 0, 0])
a = solve_lin_prog(A, b, -c)[0]
return a[-3:] if a is not None else None
Заглавная функция, принимает коэффициенты, считает случайные моменты времени, перебираем Sigma, добавляет или не добавляет выбросы, и по всему этому строит графики.
In [20]:
def solve_mnk_problem(a_real, lent=200, maxt=1):
"""
Строит графики для первой задачи.
"""
a_real = np.array(a_real)
t = np.random.uniform(0, maxt, lent)
y_real = get_y_real(a_real, t)
for pm50 in [False, True]:
for sigma in [0.1, 0.3, 0.6]:
y_corr = get_y_corr(y_real, sigma)
if pm50:
y_corr[0] += 50
y_corr[-1] -= 50
a_corr = np.empty((3, 3))
for method in range(3):
a_corr[method] = get_params(y_corr, t, method=method)
draw(a_real, a_corr, y_corr, t, maxt, sigma, pm50)
Если где-то не появилась кривая (иногда бывает с LINF), то это означает, что задача получается недопустимой (неограниченной не может быть, так как минимизируем). Это нормальная ситуация.
In [21]:
solve_mnk_problem([0.13, 0.11, 0.17], lent=50, maxt=10)
L1-норма наиболее устойчива к выбросам. При отсутствии выбросов и небольшом отклонении все нормы неплохо приближают искомую функцию.
200x80, минимизируя вектор невязок тремя разными способами.A и столбца b распределн как N(0, 1).Заглавная функция. Она же единственная и делает всю работу.
In [22]:
def solve_mnk_problem_2(m=200, n=80):
B = np.matrix(np.random.normal(0, 1, (m, n)))
c = np.matrix(np.random.normal(0, 1, (m, 1)))
E = np.eye(m, dtype=float)
e = np.ones(m, dtype=float).reshape(m, 1)
x0 = ((B.T * B).I * B.T * c).A1
A1 = np.bmat([[-E, B], [-E, -B]])
b1 = np.bmat([[c], [-c]])
c1 = np.bmat([np.ones(m), [0] * n]).A1
A2 = np.bmat([[-e, B], [-e, -B]])
c2 = np.bmat([np.ones(1), [0] * n]).A1
x1 = solve_lin_prog(A1, b1, -c1)[0]
x1 = x1[-n:] if x1 is not None else None
x2 = solve_lin_prog(A2, b1, -c2)[0]
x2 = x2[-n:] if x2 is not None else None
x0, x1, x2 = np.matrix(x0).T, np.matrix(x1).T, np.matrix(x2).T
plt.figure(figsize=(18, 6))
if x0 is not None:
plt.subplot(1, 3, 1)
plt.hist(c - B * x0)
if x1 is not None:
plt.subplot(1, 3, 2)
plt.hist(c - B * x1)
if x2 is not None:
plt.subplot(1, 3, 3)
plt.hist(c - B * x2)
plt.show()
In [23]:
solve_mnk_problem_2()
L2-норме вектор невязок распределен нормально, об этом даже говорит соответствующая теорема из курса математической статистики.L1-норме ошибки распределены по Лапласу или нормально. Не хватает теоретических знаний, чтобы сделать вывод точно.LINF-норме больше похожи на равномерное распределени. Это не очень хорошо, поскольку равномерность распределения вектора невязок означает, грубо говоря, вероятность того, что значение отклоняется на число d, такая же, как и вероятность, что оно отклоняет на d / 2.2 ** n - 1 итерацию.2^n - 1 итерацию.Как видно ниже, всего 6 переменных и совершено 63, то есть 2 ** 6 - 1, итераций.
In [24]:
A = np.matrix([[1, 0, 0, 0, 0, 0],
[4, 1, 0, 0, 0, 0],
[8, 4, 1, 0, 0, 0],
[16, 8, 4, 1, 0, 0],
[32, 16, 8, 4, 1, 0],
[64, 32, 16, 8, 4, 1]])
c = np.array([32, 16, 8, 4, 2, 1])
b = np.array([5 ** 1, 5 ** 2, 5 ** 3, 5 ** 4, 5 ** 5, 5 ** 6])
print(solve_lin_prog(A, b, c, method='lexical'))
Совершенно ясно, что и при больших n итераций также будет 2 ** n - 1. Мёожно в этом убедиться, добавляее больше переменных.
In [25]:
for n in range(7, 13):
A = np.append(A, 2 * A[-1, :], axis=0)
A[-1, -1] = 4
A = np.append(A, A[:, -1] // 2, axis=1)
A[-1, -1] = 1
c = np.append([2 ** (n - 1)], c)
b = np.append(b, [5 ** n])
print(solve_lin_prog(A, b, c, method='lexical')[2])