Лабораторная работа №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 = 20
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.05, 0.004878048780487805)

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


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)


24390.243902439026
0.0155448727568
0.0050748349777
0.0024609206183
0.00144907275689
0.000961141952353
0.000691639918162
0.000528145397028
0.000421820861725
0.000348836451365
0.000296512598602
0.000257632050157
0.000227857832477
0.000204464637565
0.000185674292204
0.000170288692657
0.000157477766924
0.00014665220181
0.000137384465095
0.000129358346155
0.000122335874094
0.000116135123416
0.000110615012639
0.000105664695807
0.000101196032101
9.7138156626e-05
9.34335097578e-05
9.00348947171e-05
8.6903270369e-05
8.40060766945e-05
8.13159509257e-05
7.88097334671e-05
7.646769107e-05
7.42729045049e-05
7.22107819585e-05
7.0268669363e-05
6.84355360846e-05
6.67017196539e-05
6.50587171013e-05
6.34990133348e-05
6.20159391597e-05
6.06035531658e-05
5.92565429496e-05
5.79701420889e-05
5.6740060023e-05
5.55624225594e-05
5.44337211756e-05
5.33507696343e-05
5.23106667082e-05
5.13107640308e-05
5.03486382669e-05
4.94220669364e-05
4.85290073413e-05
4.76675781384e-05
4.68360431747e-05
4.60327972648e-05
4.52563536429e-05
4.45053328587e-05
4.37784529282e-05
4.30745205721e-05
4.23924234049e-05
4.17311229534e-05
4.10896484037e-05
4.04670909869e-05
3.98625989283e-05
3.92753728946e-05
3.87046618809e-05
3.81497594887e-05
3.76100005514e-05
3.70847580691e-05
3.65734404194e-05
3.60754888164e-05
3.55903749903e-05
3.51175990662e-05
3.46566876221e-05
3.4207191908e-05
3.37686862106e-05
3.33407663495e-05
3.2923048293e-05
3.25151668815e-05
3.21167746493e-05
3.17275407358e-05
3.13471498782e-05
3.0975301478e-05
3.06117087357e-05
3.02560978475e-05
2.99082072585e-05
2.9567786968e-05
2.92345978834e-05
2.89084112167e-05
2.85890079226e-05
2.82761781731e-05
2.7969720867e-05
2.76694431704e-05
2.73751600875e-05
2.70866940574e-05
2.68038745767e-05
2.65265378452e-05
2.6254526433e-05
2.59876889678e-05
2.57258798409e-05
2.54689589299e-05
2.52167913384e-05
2.49692471502e-05
2.47262011973e-05
2.44875328415e-05
2.42531257676e-05
2.40228677888e-05
2.37966506615e-05
2.35743699116e-05
2.33559246687e-05
2.31412175104e-05
2.29301543132e-05
2.27226441129e-05
2.25185989708e-05
2.23179338474e-05
2.21205664824e-05
2.19264172803e-05
2.17354092026e-05
2.15474676638e-05
2.13625204344e-05
2.11804975467e-05
2.10013312061e-05
2.08249557069e-05
2.06513073511e-05
2.04803243716e-05
2.0311946859e-05
2.01461166914e-05
1.99827774677e-05
1.98218744437e-05
1.96633544711e-05
1.95071659395e-05
1.93532587202e-05
1.92015841135e-05
1.90520947974e-05
1.89047447788e-05
1.87594893472e-05
1.86162850293e-05
1.84750895468e-05
1.83358617751e-05
1.81985617039e-05
1.80631503994e-05
1.79295899683e-05
1.77978435231e-05
1.76678751482e-05
1.75396498687e-05
1.74131336191e-05
1.72882932138e-05
1.71650963191e-05
1.70435114252e-05
1.6923507821e-05
1.68050555678e-05
1.66881254756e-05
1.65726890799e-05
1.64587186187e-05
1.63461870113e-05
1.62350678372e-05
1.6125335316e-05
1.60169642885e-05
1.59099301974e-05
1.58042090698e-05
1.56997774997e-05
1.55966126314e-05
1.54946921432e-05
1.53939942319e-05
1.5294497598e-05
1.51961814311e-05
1.50990253955e-05
1.50030096176e-05
1.49081146718e-05
1.48143215689e-05
1.4721611743e-05
1.46299670404e-05
1.4539369708e-05
1.44498023821e-05
1.43612480782e-05
1.42736901802e-05
1.41871124309e-05
1.41014989221e-05
1.40168340853e-05
1.39331026827e-05
1.38502897985e-05
1.37683808304e-05
1.36873614816e-05
1.36072177523e-05
1.35279359327e-05
1.34495025953e-05
1.33719045872e-05
1.3295129024e-05
1.32191632822e-05
1.31439949931e-05
1.30696120361e-05
1.29960025329e-05
1.29231548411e-05
1.28510575486e-05
1.27796994677e-05
1.27090696301e-05
1.2639157281e-05
1.25699518743e-05
1.25014430674e-05
1.24336207165e-05
1.23664748715e-05
1.22999957719e-05
1.22341738419e-05
1.21689996862e-05
1.2104464086e-05
1.20405579944e-05
1.19772725329e-05
1.19145989873e-05
1.18525288038e-05
1.17910535855e-05
1.1730165089e-05
1.16698552205e-05
1.16101160327e-05
1.15509397215e-05
1.14923186228e-05
1.14342452092e-05
1.13767120872e-05
1.13197119943e-05
1.12632377957e-05
1.1207282482e-05
1.11518391663e-05
1.10969010813e-05
1.10424615773e-05
1.09885141188e-05
1.09350522831e-05
1.0882069757e-05
1.0829560335e-05
1.07775179167e-05
1.0725936505e-05
1.06748102036e-05
1.06241332149e-05
1.05738998381e-05
1.05241044672e-05
1.04747415889e-05
1.04258057808e-05
1.03772917097e-05
1.03291941294e-05
1.02815078791e-05
1.02342278821e-05
1.01873491433e-05
1.01408667483e-05
1.00947758613e-05
1.0049071724e-05
1.00037496535e-05
9.95880504129e-06
Out[7]:
<utils.Lab1OptCtrlModel at 0x7fdd41f79d68>

In [8]:
model.final_step


Out[8]:
245

In [9]:
X_ = np.arange(0., p_d['l'] + p_d['dh'], p_d['dh'])
Y_ = np.arange(0., p_d['T'] + p_d['dt'], p_d['dt'])

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


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.9588050412929296e-06, 0.10000995880504129)