Семинар 10

1. Материалы
2. Введение в python
3. Вариационное исчисление
4. Оптимальное управление

1. Материалы

Исходный код доступен в файлах:

intro.py - введение в python
cov.py - вариационное исчисление
oc.py - оптимальное управление

Актуальных видеозаписей нет, но если очень надо, то можно пока посмотреть старые (там мне не нравится порядок рассказа и многое из содержания, но, например, показываю, как можно делать с помощью гугл коллаба, без необходимости разворачивать среду разработки локально):

Также рекомендую посмотреть семинары Николая Петровича про решение в excel (набрать высокие баллы за ДЗ можно только по задачам на python, лучше потратить немного времени и решить максимально много задач):

2. Введение в python

2.1. Переменные, базовые типы и операции над ними

Python - язык с динамической типизацией, то есть о том, что мы передали строку в функцию, принимающую на вход целое число, узнаем только в момент исполнения кода. Справделивости ради, начиная с версии 3.5, были добавлены type hints, которые при попытки передать в функцию не то позволят редактору показать предупреждение. Но если предупреждение не было замечено по той или иной причине, то ошибка возникет только при непосредственном исполнении кода (что может быть часу на пятом-шестом его работы, к слову).

x = 30
x = True

y = 0.5
s = "Hello, Мир!"

Как видно из примера выше, переменные объявляются с помощью синтаксиса var_name = value, а тип переменной определяется значением. Обратим внимание, что в переменную можно "перезаписать" значение другого типа. Примечание со звездочкой: связано это с тем, что в python даже в случае примитивных типов переменная является указателем на объект, а поэтому ничего не "перезаписывается" в память, выделенную для значения предыдущего типа, а просто переменная начинает указывать на другое место в памяти, которое было выделено под значение нового типа. Как упражение, подумайте, что будет выведено в результате исполнения следующего кода.

a = 30
b = a

b = 42
print(a) # => 42 или 30?

Также, можем убедиться, что примитивные типы (навроде целых чисел или логического типа) обернуты в объект:

print(f"type of 1 is {type(1)}")       # => <class 'int'>
print(f"type of True is {type(True)}") # => <class 'bool'>

Далее, рассмотрим доступные нам примитивные типы и операции над ними.

2.1.1. Числовые типы

В python нам доступны числовые типы для целых чисел (int) и чисел с плавающей точкой (float).

x = 42
print(f'type of x is: {type(x)}') # => int

y = 42.0
print(f'type of y is: {type(y)}') # => float

Используем термин "число с плавающей точкой" вместо действительного числа намеренно, так как арифметика может давать неожиданные результаты в силу того, что числа с плавающей точкой хранятся в памяти компьютера с ограниченной точностью:

1.0 + 1.0 + 1.0 == 3.0 # => False

Доступны все операции над числовыми типами, которые мы можем ожидать, включая возвдение в степень и взятие по модулю:

x = 11
print(x + 31)   # => 42
print(x - 12)   # => -1
print(x * 3)    # => 33
print(x / 2)    # => 5.5
print(x // 2)   # => 5     (целочисленное деление)
print(x ** 2)   # => 121   (возведение в степень)
print(x % 3)    # => 2     (остаток от деления)

Обратим внимание, что если один из операндов float, то операнд типа int будет автоматически приведен к типу float:

print(x + 2.0)  # => 13.0
print(x ** 2.0) # => 121.0

2.1.2. Сокращенные операторы

В python в отличие от большинства С-подобных языков отсутствует оператор инкремента. Присутствуют сокращенные операторы +=, -=, *= и /=:

i = 1
while i <= 10:
    print(f'square of {i} is {i ** 2}')
    i += 1

Сокращенный оператор (на примере lhs += rhs) это сокращенный аналог (но не просто краткая запись) выражения lhs = lhs + rhs. Выражение слева (lhs) может состоять не только из имени переменной, и в случае сокращенного оператора оно вычисляется только один раз, что важно, например, если это вычисление дорогое, например, выражение слева это получение элемента списка по индексу, расчтанному с помощью дорогой функции. Также не стоит забывать о возможных ошибках при вычислении lhs дважды, если присутсвует стохастика. Ниже пример, показывающий различия:

def id_with_print(x):
    print(f"function called with x = {x}")
    return x

some_list = [0, 1, 2, 3, 4, 5]

# with +=
some_list[id_with_print(2)] += 4
# => "function called with x = 2" печатается один раз

# without +=
some_list[id_with_print(2)] = some_list[id_with_print(2)] + 4
# => "function called with x = 2" печатается дважды

2.1.3. Логический тип

В python присутствует логический тип bool, который может принимать два значения: истина (True) и ложь (False). Обратим внимание на первые заглавные буквы, ни true/false, ни TRUE/FALSE не являются валидными идентификаторами.

Операторы также вполне ожидаемые:

t = 42 > -90  # => True
f = 0 != 0    # => False

print(f"True AND False is {t and f}")  # => False
print(f"True OR  False is {t or f}")   # => True
print(f"     NOT False is {not f}")    # => True

Полезным свойством операторов and и or является вычисление по короткой схеме (не уверен, что такой перевод устоявшийся. В английском это звучит как short-circuit evaluation, либо иногда McCarthy evaluation в честь ключевого создателя семейства языков Lisp, где такие операторы, судя по всему, и появились впервые). Это значит, что для and второй операнд будет вычислен, только если первый истинен (а для or, наоборот, если ложный). Пример:

x = 0
y = 5
if x != 0 and y / x > 1:
    print(f"y is {y}")

Если в коде выше заменим and на побитовое И (оператор &), которое не обладает вышеуказанным свойством (то есть вычисляет оба операнда всегда), то аналогичный код будет падать с ошибкой ZeroDivisionError:

if x != 0 & y / x > 1:
    print(f"y is {y}")

# => ZeroDivisionError: division by zero

2.1.4. Выражение присвоения (*)

В самом начале мы вспомнили/узнали, как выглядит оператор присвоения и то, что называть присвоением происходящее не совсем точно. В духе C мы можем присвоить значение (далее как и ранее имеем в виду, что это будут указатели на один и тот же объект) сразу нескольким переменным:

x = y = z = 5

Все также в духе C можем предположить, что оператор присвоения ассоциативен справа и при выполении возвращает присвоенное значение, то есть имеем последовательность (x = (y = (z = 5)) => (x = (y = 5)) => (x = 5) => 5 и последнее значение просто не используется. Убедиться, что наша догадка неверна просто (код ниже не выполнится из-за ошибки SyntaxError: invalid syntax):

if (x = 5):
    print(x)

Но, начиная с версии 3.8 доступно выражение присвоения (assignment expression) :=, которое ведет себя почти так как мы ожидаем:

if (x := 5):
    print(x)    # => 5

Однако выше написано почти из-за отсутсвия ассоциативности справа. Чтобы первый пример с множественным присваиванием заработал нужны скобки, иначе SyntaxError:

(x := (y := (z := 5)))

Пример использования приведен в следующем разделе, а также можно обратиться к далеко не всегда хорошим паттернам использования этого свойства при разаработке на C.

2.1.5. Неявное приведение к логическому типу (*)

В python значения, не являющиеся значениями типа bool, неявно приводятся к логическому типу, если их пытаться использовать в местах, где он требуется (например как операнды для логических операторов). Правила приведения следующие:

0, 0.0, 0j, Νone, пустые списки ( [] ), кортежи ( () ), словари ( {} ), строки ( "" ) и последовательности ( range(0) ) неявно приводятся к False. Такие значения называются falsy
Все остальное (ненулевые числа, непустые строки и коллекции) неявно приводится к True. Такие значения называются truthy

Пример использования, например, такой:

# do_something() - функция, принимающая на вход подключение к базе данных
# conn - подключение к базе данных, которое при разрыве становится None
# reconnect() - функция, возвращающая восстановленное подключение к базе данных

do_something(conn or (conn := reconect()))

2.2. Строки

Отнесение типа string к примитивным некорректно, как и в принципе термин "примитивный" тип по отношению к python, ведь все равно даже в числовых и логических типах мы работаем с упакованными в объекты сущностями. Но и относить его к коллекциям (что было бы правильно, ибо за кадром мы все равно имеем массив UTF-8 симоволов, а, на самом деле, для оптимизации там существенно более сложная конструкция с глобальным словарем, хранящим ровно одну копию эквивалентной строки. Подробнее, например, здесь) мне совсем не хочется, так как строка не особо воспринимается как коллекция (на высоком уровне, а не уровне реализации).

Строки могут быть как в одинарных, так и в двойных ковычка. Различий, как в каком-нибудь php, где интерполяция происходит только внутри двойных ковычек, тут нет, как нет и символьного типа. Интерполяция (возможность внутри строки писать в фигурных ковычках любые валидные выражения, которые будут вычислены) осуществляется с помощью префикса f перед открывающей ковычкой.

some_string = "Καλημέρα ντουνιά!"
another_string = 'Привет, мир!'

x = 5
interpolated_string = f"x is {x} and x^2 is {x ** 2}"
print(interpolated_string)  # => x is 5 and x^2 is 25

Ответ на вопрос, какие ковычки использовать, звучит так: любые, главное везде одинаковые в рамках хотя бы одного репозитория. С настороженностью советовал бы относиться к советам использовать, например, одинарные ковычки, если в строке встречаются двойные. Тут достаточно вспомнить про \' и \", а не городить огород из разных ковычек:

print("I can use double quotes inside \"that way\"")
print('And \'single quotes\' also')

Для строк доступно большое количество методов, вот некоторые из них:

s = "  some strInG  "
print(s.strip())               # => "some strInG"
print(s.strip().upper())       # => "SOME STRING"

s = "ANOTHER STRING"
print(s.lower())               # => "another string"
print(s.replace("ING", "ONG")) # => "ANOTHER STRONG"
print(s.split(" "))            # => ["ANOTHER", "STRING"]

Можем обратить внимание, что вызов методов не меняет исходную строку, а возвращает результат как новую, что позволяет выстраивать "цепочки" вызовов методов, как в s.strip().upper(). В самом деле, строки в python, как в большинстве языков иммутабельны (неизменяемы).

2.3. Управляющие конструкции

2.3.1. Ветвление

Конструкция ветвления if в python выглядит стандартно, но вокруг условия круглые скобки не требуются

mark = 6.3

if mark < 4:
    print("неудовлетворительно")
elif mark < 6:
    print("удовлетворительно")
elif mark < 8:
    print("хорошо")
else:
    print("отлично")

Обращаю внимание, что if в примере сверху это statement, а не expression, то есть не возвращает никакого значения, чего иногда хотелось бы. Но if как выражение в языке тоже присутствует. Ниже пример, как его использовать для защиты от деления на ноль:

n_groups = 3
total_amount = 9000

amount_per_group = total_amount / n_groups if n_groups > 0 else total_amount

Конструкция switch, присутствующая в подовляющем большинстве C-подобных языков, здесь отсутствует, но (начиная с версии 3.10) доступен оператор match, который не только покрывает весь функционал switch, но и позволяет производить структурный pattern matching. Для интересующихся немного об этом в соотвествующем подразделе со звездочкой, а ниже пример использования в роли switch:

read_input = "q"

match read_input:
    case "q" | "quit":
        print("quiting...")
    case "r" | "run":
        print("running...")
    case cmd:
        print(f"you entered: {cmd}")

2.3.2. Циклы

В python нам доступны два вида циклов: while и for...in:

x = 10

while x > 0:
    print(f"cube of {x} is {x ** 3}")
    x -= 1

some_strings = ["foo", "bar", "baz"]
for s in some_strings:
    print(s)

Сказать особо про них нечего, кроме того, что C-подобного синтаксиса для for (например, for int i = 0; i < len(some_strings); i++ {...}) нет. Один важный момент только, что при обходе списка циклом for...in мы обходим по копии и, если нужно менять элементы списка, то делается это с помощью индекса:

xs = [0, -2, 3, 4, -1, 6, -12, -2, 5]

# напишем цикл, зануляющий отрицательные числа
for elem, idx in enumerate(xs):
    if elem < 0:
        xs[idx] = 0

print(f"non-negative xs is {xs}")

Если нужно пропустить итерацию, то используем инструкцию continue. Если выйти из цикла заранее, то break:

for num in range(100):
    if num > 12:
        print("We are not interested in numbers greater than 12")
        break

    if num % 2 == 0:
        print(f"Found an even number {num}")
        continue

    print(f"Found an odd number {num}")

2.3.3. Structural pattern matching (*)

По-русски существует термин сопоставление с образцом, но вживую я его не слышал ни разу, поэтому будем использовать только английский вариант. В python пришло это из функциональных языков программирования. Мэтчить можно списки, реализуем функцию map (которая применяет переданную функцию к каждому элементу списка и возвращает новый список):

def my_map(func, xs):
    match xs:
        case []:
            return []
        case [x, *xs]:
            return [func(x)] + my_map(func, xs)

print(my_map(lambda x: x ** 2, [0, 1, 2, 3, 4, 5]))

Можно мэтчить словари (по их отдельным полям), например:

student = {
    "programme": "Economics",
    "grade": "B",
    "name": "Ivan Petrov"
}

match student:
    case {"grade": "A", "name": name}:
        print(f"Student {name} passed exam")
    case {"grade": "B" | "C", "programme": "Economics", "name": name}:
        print(f"Student {name} passed exam (B and C allowed for Economics)")
    case {"name": name}:
        print(f"Student {name} did not pass exam")

Обращаю внимание, что не требуется, чтобы ветви покрывали все возможные случаи (как в обоих случая выше), так как match это statement и ничего не возвращает, что является существенным отличием от паттерн мэтчинга в функциональных языках. Если нужно обработать случай "все прочее", и не важно что именно прочее, используем конструкцию case _:. Например:

response = {
    "status": 200,
    "body": "Hello, World"
}

match response:
    case {"status": 200, "body": body}:
        print(f"Got response: {body}")
    case _:
        print("Request failed")

И кортежи тоже можно (заодно это пример, где ветви не покрывают все возможные случаи):

def connect_to_database():
    # do some work
    return ("ok", connection)

active_connection = None
match connect_to_database():
    case ("ok", conn):
        active_connection = conn

Более подробно можно прочесть в PEP-636.

2.4. Функции

Функции задаются с помощью ключевого слова def, если функция возвращает значение, то делаем это явно с помощью return (функция может и не возвращать ничего, к слову). Пример:

def abs(x):
    return x if x >= 0 else -x

В самом начале вспоминали про type hints, выглядит следующим образом:

def abs(x: int) -> int:
    return x if x >= 0 else -x

# модуль может быть как от целого, так и с плавающей точкой
# поэтому правильнее будет так
def abs(x: int | float) -> int | float:
    return x if x >= 0 else -x

Но повторюсь, это про удобство и документацию функций, которые мы определяем, а не защита от ошибки на уровни системы типов. Код, где мы вызываем abs("12") будет выполнен до этой строчки при запуске.

В python функции являются first-class citizens. Формально строгого определения, как и перевода на русский, для этого определения не существует, но имеется в виду, что функции можно передавать как аргумент в другую функцию, возвращать из функции и присвоить переменной. Функции, принимающие как аргумент другие функции, называют функциями старшего порядка. Также присутствует и синтаксис для анонимных функций, которые могут быть полезны именно как аргумент:

from functools import reduce

square_function = lambda x: x ** 3

sum_of_squares = reduce(
    lambda acc, el: acc + el,         # суммируем накопленное и элемент
    map(lambda x: x ** 2, range(11)), # возводим в квадрат числа от 0 до 10
    0                                 # начальное значение накопленного
)
# более идиоматический для python способ найти сумму квадратов
# sum_of_squares = sum([x ** 2 for x in range(11)])

print(f"sum of squares from 1 to 10 is {sum_of_squares}")

2.5. Коллекции

2.5.1. Списки

Список можно задать перечислением элементов, пустой список - []:

xs = [1, 3, "foo", "bar", False, [3, 4, 5]]
empty_list = []

На примере выше видно, что список может содержать элементы разных типов. В список можно добавлять элементы методом append, который изменяет исходный список. Также доступно индексирование, включая индексирование с конца с помощью отрицательных индексов, и взятие подсписка:

print(xs[0])    # => 1
print(xs[-1])   # => [3, 4, 5]
print(xs[1:3])  # => [3, "foo"]
print(xs[:2])   # => [1, 3]

2.5.2. List comprehension

На русском нет используемого перевода. Конструкция, позволяющая создавать списки на основе списков, применяя к элементам исходного списка функцию (функция обычно задается по месту):

numbers = [1, 2, 3, 4, 5, 6]
squared_numbers = [x ** 2 for x in numbers]
print(squared_numbers)  # => [1, 4, 9, 16, 25, 36]

Также есть возможность фильтровать элементы исходного списка:

squared_even_numbers = [x ** 2 for x in numbers if x % 2 == 0]
print(squared_even_numbers)  # => [4, 16, 36]

Аналогичные конструкции доступны для словарей и множеств (о них ниже), но нам они не потребуются. Для интересующихся:

# Для словарей
squares_cache = {x:(x ** 2) for x in range 10}
print(squares_cache)   # => {1: 1, 2: 4, 3: 9, 4: 16, ...}

# Для множеств
unique_squares = {x ** 2 for x in [1, 2, 2, 1, 3, 3, 2, 1]}
print(unique_squares)  # => {1, 4, 9}

2.5.3. Кортежи

Кортеж (tuple) состоит из некоторого числа значений, отделенных запятой:

request_result = ("Ok", 200, "{\"message\": \"hello\"}")

Элементы кортежа можно получать с помощью индекса, но более удобный способ это с помощью tuple destructuring (такой pattern matching на минималках), опуская неиспользуемые элементы с помощью специального идентификатора _. Например:

print(request_result[1])  # => 200

# destructuring
status, code, _ = request_result
print(f"{code} {status}") # => 200 Ok

Традиционное использование кортежей - это возврат нескольких значений из функции, например:

def summary(xs):
    # (!) ,будет отсортирован исходный список
    xs.sort()
    return xs[0], xs[-1], sum(xs) / len(xs)

some_list = [5, 3, 2, 4, 1]
min_val, max_val, avg_val = summary(some_list)
print(f"min: {min_val}, max: {max_val}, avg: {avg_val}")
# => min: 1, max: 5, avg: 3.0

2.5.4. Множества и словари

Множество (set) это коллекция, которая не может содержать дубликаты:

drinks = {"beer", "cola", "limonade"}
drinks.add("beer")
drinks.add("juice")
print(drinks)   # => {"beer", "cola", "limonade"}

Основное преимущество множеств - это дешевая (за \(O(1)\), так как они релизованы с помощью hashtable, против \(O(n)\) у списка) проверка принадлежности элемента множеству:

print("cola" in drinks)   # => True
print("coffee" in drinks) # => False

Словари (dictionary) это пары ключ-значение (в других языках встречаются термины map, hashmap, associative array, что про одно и то же, хотя почти везде нюансы с имплементациями):

student = {
    "name": "Ivan",
    "surname": "Ivanov",
    "age": 21,
    "grades": {
        "calculus": 67,
        "algebra": 58
    }
}

print(student["name"])  # => Ivan
print(student["grades"]["algebra"]) # => 58

# student["height"]  => KeyError: 'height'
print(student.get("height"))         # => None
print(student.get("height", -1))    # => -1

student["height"] = 178
student["name"] = "Alexander"

3. Вариационное исчисление

3.1. Задача 1

Численно найти допустимую экстремаль следующего функционала: \[\int_0^2 y^2 + y'^2 dt, \;\;\; y(0) = 0, \; y(2) = 1\]

3.1.1. Решение

Импортируем необходимые для решения библиотеки.

# Задача 1

from scipy.optimize import minimize
import matplotlib.pyplot as plt
import numpy as np

Задаем дискретную шкалу времени (t). Число делений выбирается исключительно эмпирическим путем: если выбрать слишком мало, то решение будет неточным, а если слишком много - будет долго считаться (и scipy вдобавок может вообще перстать его находить). Также рассчитаем шаг временной шкалы (dt), так он нам потребуется для расчета дискретного аналога производной.

num_of_fractions = 51
t = np.linspace(0, 2, num_of_fractions)
dt = t[1] - t[0]

print(f"t:      [{t[0]}, {t[1]}, {t[2]}, ..., {t[-1]}]")
print(f"dt:     {dt}")

Задаем оптимизируемый функционал в дискретном времени (f(y)). Мы используем назад смотрящую версию дискретной производной.

def f(y):
    return np.sum(y[1:] ** 2 + ((y[1:] - y[:-1]) / dt) ** 2)

Для оптимизационного алгоритма необходимо задать стартовые значения y0 (для y в каждый момент времени), а также ограничения для начального и конечного момента времени. Во все остальные моменты времени y не ограничен, поэтому описываем его двойкой (None, None). Кроме того, напоминаю о проблемах с арифметикой чисел с плавающей точко, поэтому ограничения на y(0) и y(2) накладываем не точные, а с небольшим зазором (10e-6 достаточно традиционный для таких задач).

y0 = [0.2 for _ in t]
print(f"y0:     [{y0[0]}, ..., {y0[-1]}]")

EPS = 10e-6
bounds = [(None, None) for _ in t]
bounds[0]  = (0.0 - EPS, 0.0 + EPS)
bounds[-1] = (1.0 - EPS, 1.0 + EPS)
print(f"bounds: [{bounds[0]}, {bounds[1]}, ..., {bounds[-2]}, {bounds[-1]}]")

Находим решение с помощью алогоритма L-BFGS-B (модификация оптимизационного алоритма второго порядка BFGS для задач с ограничениями, подробнее для заинтересовавшихся, например, хоть и здесь). Если решение не находится, то можно попробовать поменять стартовые значения и/или попробовать решить как задачу на максимум. Более того, если в условии написано, исследовать на экстремум, то необходимо проверить задачу и на минимум, и на максимум.

res = minimize(f, y0, method='l-bfgs-b', bounds=bounds)

В данной задаче мы знаем аналитическое решение: \[ y(t) = \frac{e^t - e^{-t}}{e^2 - e^{-2}} \]

Убедимся, что найденное численное решение с ним совпадает (построив соответствующий график).

def f_analytical(t):
  return (np.exp(t) - np.exp(-t)) / (np.exp(2) - np.exp(-2))

plt.plot(t, res.x, 'r', t, f_analytical(t), '--b')
plt.xlabel('t')
plt.ylabel('y')
plt.legend(('y численный', 'y аналитический'))
plt.savefig("./img/cov-1.png")
plt.clf()

На графике ниже видим, что решения совпадают.

3.2. Задача 2

Численно найти допустимую экстремаль следующего функционала: \[\int_0^1 y^2 + t^2 \cdot y'^2 dt \;\;\; y(0) = 1, \; y(1) = 2\]

3.2.1. Решение

Аналогично задаче один (только не знаем аналитическое решение, поэтому на графике одна линия - численная).

# Задача 2
num_of_fractions = 41
t = np.linspace(0, 1, num_of_fractions)
dt = t[1] - t[0]

def f(y):
  return np.sum(y[1:] ** 2 + (t[1:] ** 2) * (((y[1:] - y[:-1]) / dt) ** 2))

y0 = [1.5 for _ in t]
bounds = [(None, None) for _ in t]
bounds[0]  = (1.0 - EPS, 1.0 + EPS)
bounds[-1] = (2.0 - EPS, 2.0 + EPS)

res = minimize(f, y0, method='l-bfgs-b', bounds=bounds)

plt.plot(t, res.x)
plt.xlabel('t')
plt.ylabel('y')
plt.savefig("./img/cov-2.png")
plt.clf()

Получаем решение следующего вида:

3.3. Задача 3

Численно найти допустимую экстремаль следующего функционала: \[\int_0^1 y^2 + t^2 \cdot y'^2 dt \;\;\; y(0) = 1\]

3.3.1. Решение

Это задача 2, но без правой границы (просто её не задаем).

# Задача 3
num_of_fractions = 41
t = np.linspace(0, 1, num_of_fractions)
dt = t[1] - t[0]

def f(y):
  return np.sum(y[1:] ** 2 + (t[1:] ** 2) * (((y[1:] - y[:-1]) / dt) ** 2))

y0 = [1.5 for _ in t]
bounds = [(None, None) for _ in t]
bounds[0]  = (1.0 - EPS, 1.0 + EPS)

res = minimize(f, y0, method='l-bfgs-b', bounds=bounds)

plt.plot(t, res.x)
plt.xlabel('t')
plt.ylabel('y')
plt.savefig("./img/cov-3.png")
plt.clf()

Получаем решение следующего вида:

3.4. Задача 4

Численно найти допустимую экстремаль следующего функционала: \[ \int_0^1 y \cdot y'^2 dt \;\;\; y(0) = 1, \; y(1) = 4 \]

3.4.1. Решение

Аналогично предыдущим задачам. Аналитическое решение: \[ y = (7 t + 1)^{2/3} \]

num_of_fractions = 41
t = np.linspace(0, 1, num_of_fractions)
dt = t[1] - t[0]


def f(y):
    return np.sum(y[1:] * (((y[1:] - y[:-1]) / dt) ** 2), axis=0)

y0 = [0.5 for _ in t]
bounds = [(None, None) for _ in t]
bounds[0], bounds[-1] = (1.0 - EPS, 1.0 + EPS), (4.0 - EPS, 4.0 + EPS)

res = minimize(f, y0, method='l-bfgs-b', bounds=bounds)

def f_analytical(t):
    return (7 * t + 1) ** (2 / 3) 

plt.plot(t, res.x, 'r', t, f_analytical(t), '--b')
plt.xlabel('t')
plt.ylabel('y')
plt.legend(("y численный", "y аналитический"))
plt.savefig("./img/cov-4.png")
plt.clf()

Сравним численное решение с аналитическим:

3.5. Задача 5

Численно найти допустимую экстремаль следующего функционала: \[ \int_0^1 y + y''^2 dt \;\;\; y(0) = 1, \; y(1) = 4 \]

3.5.1. Решение

В этой задаче производная второго порядка, но ни на что особо это не влияет (запишем ее как производную от производной для удобства).

num_of_fractions = 41
t = np.linspace(0, 1, num_of_fractions)
dt = t[1] - t[0]

def f(y):
    dy = (y[1:] - y[:-1]) / dt
    d2y = (dy[1:] - dy[:-1]) / dt
    return np.sum(y[2:] + d2y**2)

y0 = [0.5 for _ in t]
bounds = [(None, None) for _ in t]
bounds[0], bounds[-1] = (1.0 - EPS, 1.0 + EPS), (4.0 - EPS, 4.0 + EPS)

res = minimize(f, y0, method='l-bfgs-b', bounds=bounds)

plt.plot(t, res.x, 'r')
plt.xlabel('t')
plt.ylabel('y')
plt.savefig("./img/cov-5.png")
plt.clf()

Получаем следующее решение

4. Оптимальное управление

4.1. Задача 1

Численно решить задачу оптимального управления: \[ \int_0^2 \frac{1}{2} y^2 + y^2 \cdot t^2 dt \] \[ y' = u, \, y(0) = 1, \, u \in [-1, 1] \]

4.1.1. Решение

Подключаем необходимые библиотеки:

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import os

Формулируем задачу (через фиктивную переменную \(z\) задаем подинтегральное выражение как \(z'\)):

# Задача 1

# Инициализируем модель
m = GEKKO(remote=False)
nt = 101
m.time = np.linspace(0,2,nt)

# Задаём переменные
y = m.Var(value=1)
z = m.Var(value=5)
u = m.Var(value=0,lb=-1,ub=1)
t = m.Var(value=0)
# Отмечаем последнюю точку временной шкалы
p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)


# Задаём уравнения
m.Equation(y.dt() == u)
m.Equation(t.dt() == 1)
m.Equation(z.dt() == 0.5 * y ** 2 + y**2 * t**2)

Задаем целевую функцию (значение функционала (в нашем случае интеграла) в последний момент времени) и запускаем оптимизатор:

m.Obj(z * final)    # Целевая функция
m.options.IMODE = 6 # Задача оптимального управления/динамического программирования
m.solve()

Рисуем получившееся решение (и состояние, и управление):

plt.plot(m.time[1:], y.value[1:], 'k-', label=r'$y$')
plt.plot(m.time[1:],u.value[1:],'r--',label=r'$u$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Value')
plt.savefig("./img/oc-1.png")
plt.clf()

Получили следущее решение:

4.2. Задача 2

Численно решить задачу оптимального управления: \[ \int_0^2 -y + \frac{1}{2} u^2 + u dt \] \[ y' = y + u, \, y(0) = 2, \, u \in \mathbb{R} \]

4.2.1. Решение

Аналогично предыдущей задаче:

# Задача 2

m = GEKKO(remote = False)
nt = 301
m.time = np.linspace(0,2,nt)

y = m.Var(value=2)
z = m.Var(value=0)
u = m.Var(value=0)
p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)


m.Equation(y.dt()==y + u)
m.Equation(z.dt()==-y + (u ** 2) * 0.5 + u)
m.Obj(z*final)
m.options.IMODE = 6
m.solve(disp=False)

В данной задаче мы знаем аналитическое решение: \[u(t) = e^{2-t}, \quad y(t) = 2 + \frac{e^2}{2} \left( e^t - e^{-t} \right) \]

Нарисуем решение, сравнив с аналитическим:

def u_a(t):
    return [np.exp(-elem + 2) - 2 for elem in t]

def y_a(t):
    return [2 + 0.5 * np.exp(2) * (np.exp(elem) - np.exp(-elem)) for elem in t]

plt.figure(figsize=(12, 10))
plt.subplot(2, 1, 1)
plt.plot(m.time, y.value, 'k-', label=r'$y$')
plt.plot(m.time, y_a(m.time), 'r--', label=r'$y^*$')
plt.legend(loc='best')
plt.ylabel('Value')
plt.subplot(2, 1, 2)
plt.plot(m.time[1:], u.value[1:], 'k-', label=r'$u$')
plt.plot(m.time[1:], u_a(m.time)[1:], 'r--', label=r'$u^*$')
plt.legend(loc='best')
plt.ylabel('Value')
plt.xlabel('Time')
plt.savefig("./img/oc-2.png")
plt.clf()

Получили следующее решение:

4.3. Задача 3

Численно решить задачу оптимального управления: \[ \int_0^2 -y + \frac{1}{2} u^2 + u dt \] \[ y' = y + u, \, y(0) = 2, \, u \in [e-2, +\infty) \]

4.3.1. Решение

Аналогично предыдущей задаче, изменилось только множество допустимых значений управления:

# Задача 3
m = GEKKO(remote = False)
nt = 301
m.time = np.linspace(0,2,nt)

y = m.Var(value=2)
z = m.Var(value=0)
u = m.Var(value=np.e-2, lb=np.e-2)
p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)

m.Equation(y.dt()==y + u)
m.Equation(z.dt()==-y + (u ** 2) * 0.5 + u)
m.Obj(z*final)
m.options.IMODE = 6
m.solve(disp=False)

Знаем аналитическое решение. При \(t \in [0,1]\): \[u(t) = e^{2-t} - 2, \quad y(t) = 2 + \frac{e^2}{2} \left( e^t - e^{-t} \right) \]

При \(t \in (1, 2]\): \[u(t) = e - 2, y(t) = 1 + \frac{e^2 - 1}{2} e^t - e + 2 \] \[\] Нарисуем и сравним с аналитическим решением:

def u_a(t):
    return [np.exp(-elem + 2) - 2 if elem <= 1
            else np.e - 2
            for elem in t]

def y_a(t):
    return [2 + 0.5 * np.exp(2) * (np.exp(elem) - np.exp(-elem)) if elem <= 1
            else (1 + 0.5 * (np.e ** 2 - 1)) * np.exp(elem) - np.e + 2
            for elem in t]

plt.figure(figsize=(12, 10))
plt.subplot(2, 1, 1)
plt.plot(m.time, y.value, 'k-', label=r'$y$')
plt.plot(m.time, y_a(m.time), 'r--', label=r'$y^*$')
plt.legend(loc='best')
plt.ylabel('Value')
plt.subplot(2, 1, 2)
plt.plot(m.time[1:], u.value[1:], 'k-', label=r'$u$')
plt.plot(m.time[1:], u_a(m.time)[1:], 'r--', label=r'$u^*$')
plt.legend(loc='best')
plt.ylabel('Value')
plt.xlabel('Time')
plt.savefig("./img/oc-3.png")