Лабораторная работа №1

по предмету "оптимальное управление"

Выполнил: студент гр. А-14м-16 Мигаль И.А
Проверил: Зубков Павел Валерьевич

In [ ]:
import sys
sys.path.insert(0, '/home/ivmig/OneDrive/Documents/A-14m-16/2CURSE3SEM/optimal_control/lab_1/')
print(sys.path)

In [1]:
print(__doc__)

# Author: Ivan Migal ivan.migal@mail.ru
# License: BSD 3 clause

import math
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from pylab import rcParams
from matplotlib import colors as mcolors
import utils

from utils import array
from utils import TDMA


colors = list(mcolors.BASE_COLORS)
colors += list(mcolors.CSS4_COLORS.keys())


Automatically created module for IPython interactive environment

In [2]:
%matplotlib inline
rcParams['figure.figsize'] = 12, 12
plt.style.use('ggplot')
# Настройка шрифта
font = {'family' : 'DejaVu Sans',
        'weight' : 'bold',
        'size'   : 16}

matplotlib.rc('font', **font)

Инициализация данных

    Обозначения см. в описании Л.Р. №1

In [3]:
# Словарь параметров
p_d = {}

In [4]:
# Заданные положительные величины
p_d['a'], p_d['l'], p_d['v'], p_d['T'] = 1., 1., 1., 2.

# Решение тестового примера
def x(s, t):
    return math.sin(t) + math.sin(s + math.pi / 2.)

# Плотность источников тепла
def f(s, t):
    return math.cos(t) + p_d['a'] ** 2 * math.sin(s + math.pi / 2.)
    
# Температура внешней среды
def p(t):
    #return 1. / p_d['v'] * math.cos(p_d['l'] + math.pi / 2) + math.sin(t) + math.sin(p_d['l'] + math.pi / 2)
    return 0.
    
# Распределение температуры в начальный момент времени
def fi(s):
    return math.sin(s + math.pi / 2.)

p_d['p(t)'] = p

p_d['f(s, t)'] = f

# Заданные числа
p_d['p_min'], p_d['p_max'], p_d['R'] = -100., 100., 100.

p_d['fi(s)'] = fi

# Желаемое распределение температуры
def y(s):
    #return s * math.sin(s + math.pi / 2.)
    return math.sin(p_d['T']) + math.sin(s + math.pi / 2.)

p_d['y(s)'] = y

# Число точек на пространственной и временной сетке соответственно
N = 30
p_d['N'], p_d['M'] = N, 10 + N * N

# Шаг на пространственной и временной сетке соответственно
p_d['dh'], p_d['dt'] = p_d['l'] / p_d['N'], p_d['T'] / p_d['M']
p_d['l'], p_d['T'], p_d['dh'], p_d['dt']


Out[4]:
(1.0, 2.0, 0.03333333333333333, 0.002197802197802198)

Создание модели и управление


In [5]:
model = utils.Lab1OptCtrlModel(p_d)

In [6]:
criterion = utils.criterion_3
get_alpha = utils.get_alpha_5

In [7]:
model.solve(criterion, get_alpha, eps=10**-5)


10989.010989010989
0.0159073299676
0.00518208265573
0.00250111580825
0.00146503326042
0.000966916037968
0.000692831701906
0.00052726912928
0.00042008983248
0.000346861114426
0.000294604520106
0.000255947363482
0.000226469214328
0.000203399154089
0.000184934023274
0.000169862271837
0.000157346954827
0.000146795664863
0.000137779945863
0.000129983892016
0.000123170508297
0.000117159184116
0.000111810294584
0.000107014476008
0.000102685029047
9.8752453313e-05
9.51604590219e-05
9.18630180499e-05
8.88221568523e-05
8.60062858607e-05
8.33889215831e-05
8.09476994366e-05
7.86636041111e-05
7.6520364315e-05
7.45039729007e-05
7.26023034622e-05
7.0804801779e-05
6.91022357816e-05
6.74864916156e-05
6.59504062718e-05
6.44876294168e-05
6.30925086876e-05
6.17599939572e-05
6.04855570254e-05
5.92651239213e-05
5.80950175708e-05
5.69719090264e-05
5.58927758041e-05
5.48548661453e-05
5.38556682425e-05
5.28928836389e-05
5.19644041539e-05
5.1068291797e-05
5.02027612266e-05
4.93661643805e-05
4.85569769696e-05
4.77737865732e-05
4.70152821159e-05
4.62802445412e-05
4.55675385238e-05
4.48761050853e-05
4.42049550005e-05
4.3553162894e-05
4.29198619444e-05
4.2304239122e-05
4.17055308983e-05
4.11230193716e-05
4.05560287621e-05
4.00039222358e-05
3.94660990192e-05
3.89419917762e-05
3.84310642167e-05
3.79328089148e-05
3.74467453133e-05
3.69724178983e-05
3.65093945233e-05
3.60572648721e-05
3.56156390439e-05
3.51841462517e-05
3.4762433621e-05
3.43501650815e-05
3.39470203423e-05
3.35526939439e-05
3.31668943789e-05
3.27893432775e-05
3.24197746505e-05
3.20579341855e-05
3.17035785921e-05
3.13564749924e-05
3.10164003521e-05
3.06831409494e-05
3.03564918798e-05
3.00362565919e-05
2.97222464542e-05
2.94142803479e-05
2.91121842869e-05
2.88157910596e-05
2.8524939894e-05
2.82394761421e-05
2.79592509835e-05
2.76841211466e-05
2.74139486454e-05
2.71486005327e-05
2.6887948666e-05
2.66318694875e-05
2.63802438161e-05
2.61329566511e-05
2.5889896986e-05
2.56509576331e-05
2.5416035057e-05
2.51850292174e-05
2.49578434192e-05
2.47343841718e-05
2.45145610542e-05
2.42982865878e-05
2.40854761155e-05
2.38760476865e-05
2.36699219476e-05
2.34670220383e-05
2.32672734931e-05
2.30706041462e-05
2.2876944043e-05
2.26862253539e-05
2.24983822929e-05
2.23133510408e-05
2.21310696701e-05
2.1951478075e-05
2.17745179038e-05
2.16001324942e-05
2.1428266812e-05
2.12588673923e-05
2.10918822827e-05
2.09272609901e-05
2.07649544287e-05
2.06049148709e-05
2.04470959e-05
2.0291452365e-05
2.0137940337e-05
1.99865170681e-05
1.98371409509e-05
1.96897714811e-05
1.95443692199e-05
1.94008957594e-05
1.92593136888e-05
1.91195865614e-05
1.89816788639e-05
1.88455559865e-05
1.87111841935e-05
1.85785305963e-05
1.84475631262e-05
1.8318250509e-05
1.81905622406e-05
1.80644685625e-05
1.79399404399e-05
1.78169495388e-05
1.76954682055e-05
1.75754694457e-05
1.74569269049e-05
1.73398148496e-05
1.72241081489e-05
1.71097822566e-05
1.69968131946e-05
1.68851775361e-05
1.67748523899e-05
1.66658153851e-05
1.65580446563e-05
1.64515188292e-05
1.63462170068e-05
1.62421187565e-05
1.61392040966e-05
1.60374534841e-05
1.5936847803e-05
1.5837368352e-05
1.57389968338e-05
1.5641715344e-05
1.55455063606e-05
1.54503527335e-05
1.53562376752e-05
1.52631447508e-05
1.5171057869e-05
1.50799612728e-05
1.49898395314e-05
1.4900677531e-05
1.48124604675e-05
1.4725173838e-05
1.46388034332e-05
1.45533353304e-05
1.44687558858e-05
1.43850517279e-05
1.43022097503e-05
1.42202171058e-05
1.41390611992e-05
1.4058729682e-05
1.39792104454e-05
1.39004916156e-05
1.38225615471e-05
1.37454088179e-05
1.3669022224e-05
1.35933907738e-05
1.35185036838e-05
1.34443503731e-05
1.33709204591e-05
1.32982037524e-05
1.32261902527e-05
1.31548701442e-05
1.30842337918e-05
1.30142717361e-05
1.29449746903e-05
1.28763335359e-05
1.28083393187e-05
1.27409832455e-05
1.26742566803e-05
1.26081511408e-05
1.25426582951e-05
1.24777699581e-05
1.24134780887e-05
1.23497747863e-05
1.22866522881e-05
1.22241029656e-05
1.21621193221e-05
1.21006939899e-05
1.20398197273e-05
1.1979489416e-05
1.19196960586e-05
1.18604327758e-05
1.18016928042e-05
1.17434694939e-05
1.16857563056e-05
1.16285468091e-05
1.15718346804e-05
1.15156136998e-05
1.14598777497e-05
1.14046208127e-05
1.13498369691e-05
1.12955203954e-05
1.12416653623e-05
1.11882662325e-05
1.1135317459e-05
1.10828135835e-05
1.10307492346e-05
1.09791191257e-05
1.09279180538e-05
1.08771408975e-05
1.08267826158e-05
1.07768382462e-05
1.07273029033e-05
1.06781717773e-05
1.06294401326e-05
1.05811033062e-05
1.05331567067e-05
1.04855958124e-05
1.04384161705e-05
1.03916133955e-05
1.0345183168e-05
1.02991212333e-05
1.02534234008e-05
1.02080855419e-05
1.01631035896e-05
1.01184735371e-05
1.00741914366e-05
1.00302533982e-05
9.98665558924e-06
Out[7]:
<utils.Lab1OptCtrlModel at 0x7fba4bb35e10>

In [8]:
model.final_step


Out[8]:
261

In [9]:
# Здесь arange странно
X_ = np.arange(0., p_d['l'] + p_d['dh'] / 2., p_d['dh'])
Y_ = np.arange(0., p_d['T'] + p_d['dt'], p_d['dt'])
len(X_), len(model.y_arr), p_d['dh'] * 31


Out[9]:
(31, 31, 1.0333333333333332)

Процесс управления


In [10]:
y_s = model.y_arr
bounds = [min(min(y_s), min(model.x_arr[-1][-1,:])), max(max(y_s), max(model.x_arr[-1][-1,:]))]

In [11]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(X_, y_s, color='r', label='Желаемое')
part = 1
count = 5
step = int(model.final_step / count / part)
for i in range(0, int(model.final_step / part), step):
    ax.plot(X_, model.x_arr[i][-1,:], color=colors[int(i / step)], linestyle='--', label='{}-я итерация'.format(i))
ax.plot(X_, model.x_arr[-1][-1,:], color='b', label='Последняя итерация')
ax.set_ylim(bounds)
plt.xlabel('s')
plt.ylabel('y(s)')
plt.title('Распределения температуры стержня в процессе управления')
plt.legend(title='Легенда', loc='center left', bbox_to_anchor=(1.0, 0.5))
plt.show()



Желаемое и полученное распределение температур


In [12]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(X_, y_s, color='r', label='Желаемое')
ax.plot(X_, model.x_arr[-1][-1,:], color='b', label='Полученное')
ax.set_ylim(bounds)
plt.xlabel('s')
plt.ylabel('y(s)')
plt.title('Желаемое и полученное распределения температуры стержня')
plt.legend(title='Легенда', loc='center left', bbox_to_anchor=(1.0, 0.5))
plt.show()



Управление


In [13]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(Y_, model.p_arr[0], color='b', label='Начальный момент времени')
ax.plot(Y_, model.p_arr[-1], color='r', label='Полученное после управления')
plt.xlabel('t')
plt.ylabel('p(t)')
plt.title('Управление начальный и окончательный момент')
plt.legend(title='Температуры внешней среды', loc='center left', bbox_to_anchor=(1.0, 0.5))
plt.show()



График изменения ошибки


In [14]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(model.err, color='black', label='Ошибка')
ax.set_ylim([min(model.err), min(model.err) + .1])


Out[14]:
(9.9866555892350993e-06, 0.10000998665558924)