CC BY
59
K = load('K.dat'); M = load('M.dat');
%определение числа строк в матрице K n = size(K,1);
%ввод единичной I и нулевой Z матриц размерности n
I = eye(n); Z = zeros(n); %формирование матриц A и B
A = vertcat(horzcat(Z,I),horzcat(-
inv(M)*K,Z)); B = vertcat(Z,inv(M));
%размерность уравнения Риккати dim, матрицы E и S
dim = size (A, 1); E = eye (dim) ; S = B*transpose(B);
%начальные значения P и ошибки P = zeros(dim); SE(1) = norm(CalcFR(S, A, E, P));
^допустимая ошибка minErr, шаг по времени h minErr = 0.001; h = 0.1; %максимальное число итераций maxIt = 100;
%начальный номер итерации it = 1;
while ((SE(it)>minErr)&(it<maxIt)) it = it + 1;
%расчет следующего приближения к P F1 = h*CalcFR(S, A, E, P); F2 = h*CalcFR(S, A, E, P + F1/2); F3 = h*CalcFR(S, A, E, P + F2/2); F4 = h*CalcFR(S, A, E, P + F3); P = P + (F1 + 2*F2 + 2*F3 + F4)/6; SE(it) = norm(CalcFR(S, A, E, P)); end;
Исследуем эффективность приведенного кода на примере расчета управления конечно-элементным объектом, который состоит из последовательно соединенных n 20-узловых изопараметрических элементов (рис. 1). Каждый узел считаем наделенным тремя степенями свободы (перемещения в трех направлениях декартовой системы координат).
Результаты тестирования
/ — Г ) —о- г / А
/ —£>— у -ft— -< / -в— - У -ft— -А У
Рисунок 1 - Исследуемый конечно-элементный объект
В таблице 1 приведены данные по быстродействию разработанного кода, полученные на процессоре Intel Celeron частотой 2.2 ГГЦ с ОЗУ 2 Гб под управлением Windows XP.
кода MRDE.m. Таблица 1
Рисунок 2 - Время синтеза оптимального управления
Число элементов Число узлов Размерность K и M Размерность P Время расчета P, сек
1 20 60 x 60 120 x 120 2,929936
2 32 96 x 96 192 x 192 10, 4 67 629
3 44 132 x 132 264 x 264 26,467749
4 56 168 x 168 336 x 336 52,614708
5 68 204 x 204 408 x 408 98,437217
6 80 240 x 240 480 x 480 162,432924
7 92 276 x 276 552 x 552 251,138686
8 104 312 x 312 624 x 624 357,065507
9 116 348 x 348 6 9 6 x 6 9 6 501,303477
10 128 384 x 384 768 x 768 677,600257
Результатам, приведенным в таблице 1, соответствует график на рис. 2.
Полученные данные свидетельствуют о высокой производительности разработанного методического комплекса.
ЛИТЕРАТУРА
1. Партон В. З., Кудрявцев Б. А. Электромагнитоупругость пьезоэлектрических и электропроводных тел. - М.: Наука, 1988, 472 с.
2. Горшков А. Г., Рабинский Л. Н., Тарлаковский Д. В. Основы тензорного анализа и механика сплошной среды. Учебник для вузов. - М.: Наука, 2000, 214 с.
3. О. Зенкевич. Метод конечных элементов в технике: Пер. с англ.- М.: Мир, 1975, 541 с.
4. M. Shafiee, S. Amani. Optimal control for a class of singular systems using neural network. Iranian Journal of Science & Technology, Transaction B, Engineering, Vol. 29, No. B1, 2005, 3348.
5. Горячев Н.В. Исследование и разработка средств и методик анализа и автоматизированного выбора систем охлаждения радиоэлектронной аппаратуры / Горячев Н.В., Танатов М.К., Юрков Н.К. // Надежность и качество сложных систем. 2013. № 3. С. 70-75.
6. Кочегаров И.И. Программный пакет моделирования механических параметров печатных плат / Кочегаров И.И., Таньков Г.В. // Труды международного симпозиума Надежность и качество. 2011. Т. 2. С. 334-337.
7. Сабоннадьер Ж.-К., Кулон Ж.-Л. Метод конечных элементов и САПР: Пер. с франц. - М.: Мир, 1989, 190 с.
УДК 62-9, 621.3.078 Ерофеев С.А.
ФГБОУ ВПО «Санкт-Петербургский государственный политехнический университет», Санкт-Петербург, Россия
РАСЧЕТ ДИНАМИЧЕСКИХ СИСТЕМ И УПРАВЛЕНИЕ ИМИ С ПОМОЩЬЮ НЕЙРОННЫХ СЕТЕЙ НА БАЗИСНЫХ ПРОБНЫХ ФУНКЦИЯХ В СРЕДЕ MATLAB
Динамический объект, как правило, описывает- ференциальных уравнений состояния. Изоморфизм ся дифференциальным уравнением или набором диф- [1] этих уравнений, то есть их аналогичность
для ряда механических, электрических, гидравлических и других систем, позволяет применять одни и те же методы разрешения этих уравнений с целью расчета выходных параметров. Выбор этих методов часто исходит из невозможности аналитического получения интересующих данных, поэтому достаточно распространенной концепцией является применение хорошо зарекомендовавших себя численных алгоритмов, допускающих высокоэффективную программную реализацию. К числу таких технологий относятся нейронные сети различных разновидностей, среди которых выделим и исследуем многослойный персептрон, работающий на базисных пробных функциях [2]. Для работы с сетями такого типа необходимы адаптивное алгоритмическое описание процедуры обучения сети и технология программной реализации алгоритма. Восполним фактический недостаток такого рода информации.
Рассмотрим объект, поведение которого описывается дифференциальным уравнением первого порядка
^ = Г (1,у)
(1)
1 п
Е=2 ^
2 к=1
¿у
(к)
¿г
--Г
(к)
где у(к) = у(гк ) , Г(к) = Г (гк ) , должна стать меньше определенного значения minErr, то есть должна быть минимизирована. Значения гк будут распределяться равномерно, исходя из их общего числа п на интервале [а, Ь]. Нейронная сеть будет итерационным образом уточнять веса и смещения по правилу обучения [3]:
Ш(к) = Ш(к) дЕ В(к) = В(к) дЕ «1 п д , Bj Bj п
3В« ,
1 = V® -п
дЕ д1 .
(3)
где
№(В , V® , В( 1 , 1 , 1
(к)
веса и смещение нейрона j
коэффициент
с неизвестной величиной у скалярного типа. Такое уравнение, в частности, может отражать механическое поступательное или вращательное движение в одном направлении, заполнение резервуара гидравлической жидкостью или ее отвод, растворение вещества и другие изоморфные процессы, а под у соответственно может пониматься линейная или угловая скорость тела, уровень жидкости, концентрация растворенного вещества и т. д. [1] .
Решение задачи (1) с начальным условием УОо) = Уо в соответствии с методом пробных функций находится в виде у(1) = Уо + (г — го)N(0 , где N(t) - выходное значение нейронной сети типа трехслойного персептрона:
сети, в которую подается гк ; п обучения.
Примем г0 = 0 , у0 = 0 , что влечет у(1) = гЫ(г) .
Зададимся униполярной сигмовидной функцией ак-
1
тивации о(х) = - — 1 + е 1
вых частях выражений (1) Вначале рассмотрим
и найдем производные в пра-
дЕ дW1®
с учетом (2):
дЕ д W®
VI — гт 1! д (д У« 1 а Г (,) | дг, 1 ДаW1(k) [ дг, ) дwf)j
(4)
Так как у = г^^) и ^(х) = о(х)(1 — а(х)) , то ¿х
ду<" -- N(1,) + г,= 5 V® 2® + г, V V® 2®(1—г®уw® (5)
дг,
дг,
1=1
1=1
где 7((,) = о(г^'' + В®)
Выполнив ряд преобразований, получим: д ( ду(° ^
дW1(k) 1 дг(
= 1кк/к) (1 — г )к) )[1к^(к) (1 — 2 2 |к)) + 2]
7 (к)
(к)(
7(к))
где 7((к) =*(гк1 + В®) .
Затем преобразуем:
д/(,) _ дГ® дут
(к)
= = г2¥® 7(к)С1 — 7(к)) дГ'
дW1(k) ду(0 дW1(k) к ^ 71 (1 71 ) ду(к)
(6)
(7)
ш --- )
Рисунок 1 - Трехслойный персептрон
где наряду с одним входным сигналом t на входном слое и одним выходным сигналом N на выходном слое присутствует скрытый слой нейронов 1 ... т. Каждый нейрон j выполняет роль преобразователя входного сигнала, умножая его на вес Wj ,
добавляя смещение В^ и применяя сигмовидную функцию активации о, что приводит к получению выходного сигнала нейрона Zj = (т(гWj + Bj) . Эти
сигналы умножаются на веса Vj , складываются в
выходном слое, в результате чего формируется т
выходное значение N = 5 VjCT(fWj■ + Bj) .
1=1
Поставим задачу подбора таких значений весов ( Wj , Vj ) и смещений В^ , при которых соотношение (1) будет выполняться с указанной точностью на заданном интервале 1 е [а,Ь] . Для этого в соответствии с концепцией метода пробных функций суммарная ошибка
Таким же образом можно вывести следующие соотношения:
д 1 = V® 2® (1 — 2® )[г^®(1 — 22®) + 1],
дВ(к) I дг,
дГ
(к) = г12,®(1 — 2(к))- ... дВ(к) к 1 1 1 д(к)
7(к)) д ду1
(к)
(8)
дУ® I дг,
^ 1 = 2®^®(1 — 21®) + 1],
дГ('
- = г2(к)
дГ1
(к)
дУ® к 1 ду(к)
1
(9)
дт(о дГ(°
Заметим, что -;л- = 1к—тт-
, дW1(k) к дВ((к)
Производные
дЕ дЕ
-^ и —^ по форме аналогичны
дV(k) дВ((к) ф р
(4), а полученные выражения (8) и (9) позволяют однозначно рассчитывать все правые части уравнений (3).
Рассмотрим задачу (с г0 = 0 и у0 = 0 ) :
ду = 2 у — у2 +1 ,
а ,
(10)
которая имеет точное аналитическое решение [4]:
y = 1 + л/2 tgh
S t+1
2 I V2 +1
Решим это уравнение с помощью представленной выше нейронной сети, выполнив ее программную реализацию на языке MATLAB. Вначале выделим функцию TPoint, позволяющую на интервале [a, b] по индексу k и числу точек n определить соответствующий дискретный момент времени tk :
function [t] = TPoint(k,n) a = 0; % начало временного интервала b = 1; % конец временного интервала t=(b-a)/(n-1)*(k-1); end
Кроме того, оформим вычисление epsk = — f
dt
и yk в виде отдельной функции:
function [eps,Y] = EpsY(W,V,B,t,m) N = 0; dN = 0; for j=1:m
arg = W(j)*t + B (j) ; Z = 1/(1 + exp(-
arg));
N = N + V (j ) * Z; dN = dN + Z*(1-Z)*V(j)*W(j); end;
dY = N + t*dN; Y = t*N; Fr = 2*Y - YA2 + 1; eps = dY - Fr; end
Обновление весов согласно (3)-(9) и вычисление суммарной ошибки (2) также разместим в отдельной функции (learnRate - коэффициент обучения (= п)):
function [err, nW,nV,nB] =
NetTrain(n,learnRate,W,V,B,tn,m)
newW = W; newB = B; newV = V; for j = 1:m
dEdV = 0; dEdW = 0; dEdB = 0; for i=1:n t=tn(i);
[eps,Y] = EpsY(W,V,B,t,m); arg = t*W(j) + B (j); Z = 1/(1+exp(-Z*(1-Z);
arg)); dZ
2*Z)+2);
ddYdB ddYdW
ddYdV dFrdY
t*V(j)*dZ*dFrdY;
dFrdW
V(j)*dZ*(t*W(j)*(1-2*Z)+1); = t*V(j ) *dZ* (t*W (j ) * (1-
Z*(t*W(j)*(1-Z)+1); = 2 - 2*Y; dFrdB =
t*dFrdB; dFrdV = t*Z*dFrdY dEdB + eps*(ddYdB - dFrdB) dEdW + eps*(ddYdW - dFrdW) dEdV + eps*(ddYdV - dFrdV)
dEdB dEdW dEdV end;
%рассчитываем новые значения весов
V(j) - learnRate*dEdV; W(j) - learnRate*dEdW; B(j) - learnRate*dEdB;
newB;
смещения
newV(j) newW(j) newB(j) end;
nW = newW; nV = newV; nB %рассчитываем ошибку err = 0; for i = 1:n
t = TPoint(i,n); [eps, Y] = EpsY(nW,nV,nB,t,m); err = err + eps*eps/2; end;
end
С учетом введенных функций главный модуль нейронной сети выглядит так:
maxIt = 5000; %максимальное количество итераций
n = 6; %число дискретных моментов времени m = 6; %число нейронов в скрытом слое minErr = 0.001; %минимальная ошибка learnRate = 0.05; %коэффициент обучения %инициализация весов и смещений W = rand (m,1) - 0.5; B = rand(m,1) - 0.5; V = rand(m,1) - 0.5;
%рассчитываем начальную ошибку err = 0;
for i = 1:n
tn (i) = TPoint(i,n); [eps,Y] =
EpsY(W,V,B,tn(i),m); err = err + eps*eps/2; end; it = 0;
%в цикле обновляем веса и смещения, вычисляем суммарную ошибку err
while ((err>minErr)&(it<maxIt))
it = it + 1; [err,nW,nV,nB] = NetTrain(n,learnRate,W,V,B,tn,m); disp(err); W = nW; V = nV; B = nB; end;
Условием выхода из цикла, то есть отмены повторного прогона нейронной сети, является либо получение ошибки меньше минимального значения miner, либо превышение максимального количества итераций maxIt (во избежание зацикливания сети). Начальные значения весов и смещений инициализируются случайным образом с учетом правила Нгуена-Видроу [3]. Оператор disp выводит ошибку на экран с целью контроля сходимости.
Приведенный код был протестирован на процессоре Intel Celeron частотой 2.2 ГГЦ с ОЗУ 2 Гб под управлением Windows XP и показал частую сходимость в пределах 1-2 секунд, что является признаком высокой производительности по сравнению с альтернативными методами, например, методом Рунге-Кутта. С помощью встроенных средств визуализации MATLAB выполнено сравнение теоретических и нейросетевых расчетных данных в графическом виде (рис. 1), которое доказывает правильность получаемы:': результатов.
Рисунок 2 - Теоретические (-) и нейросетевые (о) результаты решения задачи (10)
Теперь рассмотрим перспективы разработанной технологии применительно к задаче управления. Известно [5], что поиск оптимального управления, например, в рамках линейно-квадратичного регулятора может выполняться посредством решения уравнения Риккати, которое в несингулярном
дР
случае сводится [6] к — = £(Р) , где Р - матри-
дг
ца. Такое матричное уравнение эквивалентно системе обычных скалярных дифференциальных уравне-дРг
ний вида —- = А(/,р,) , где 1 и j дt
j'
номер строки
и номер столбца в матрице Р. В этой связи распространим технологию реализованной выше нейронной сети на решение системы дифференциальных уравнений:
^=А(и у„ У2>-> Ур)
dy2_ dt
--f2(t, Л У2,-, Ур )
dyp и
"dt1 = fp(t- У1, У2>-> Ур)
, (11)
Для этого достаточно обобщить ошибку (2) до
dt
E =
1 n p
1Е Е
2
k=1 q=1
f (k) Jq
(12)
и вместо одной нейронной сети создать р трехслойных персептронов (рис. 1). Выходные значения персептронов согласно (12) будут складываться, что по концепции (3) приводит к взаимосвязи нейронных сетей по весам и смещениям.
и
В качестве примера рассмотрим систему ураЕ нений:
которая
точное
dyi - у dt У2
^ - 1 - у, dt л
решение
(13)
Yi(i) -1 - cos(t)
у 2 (г) = ^"(г) .
Обновление весов и смещений будет выполняться по соотношениям, аналогичным (3), то есть
W1(k) - W1(k)-n-^ , B1(k) - B1(k)-n J J /5W1(k) , 1 1
V1(k) - V1(k) -n
B2
(k) .
B2(k
dE 5V1f '
dE
W2((k) - W2(k) - n
(k).
dE 5B1((k)
dE
dW2
(k)
1 ""J '' SB2f
V2(k) - V2(k) -n 8E
5V2((k) '
где 111, VI, В1 - веса и смещения первого пер-септрона, а 112, V2, В2 - веса и смещения второго персептрона. Обозначив правые части первого и второго уравнений как ^ и :2 , выразим произ-
(к) .
водную, входящая в выражение для W1( dE
3W1®
-№ f д | ду1 f + iM f d (д f -il dt - Iaw1®I dt ~f1 I dt ~f2 Iaw1®I dt ~f2
-11
учетом
f2 - 1-У
понятно,
д U
f? I- дf2
Ф0 , то есть веса первой
дWlf) ^дг 72 J дWlf)
и второй сетей взаимосвязаны. Дальнейшие преобразования аналогичны выполненным выше выкладкам для уравнения (1). Появление дополнительных составляющих и перекрестных взаимосвязей приводит к расширению программной реализации. В частности, производные функции ошибки (12) по весам и смещениям вычисляются теперь следующим образом:
dEdW1=dEdW1+eps(1)*ddY1dW1 + ерз(2^УЫ1Т1; dEdW2=dEdW2+eps(2)*ddY2dW2 - eps(1)*dY2dW2; dEdB1=dEdB1+ eps(1)*ddY1dB1 + eps(2)*dY1dB1 dEdB2=dEdB2+ eps(2)*ddY2dB2 - eps(1)*dY2dB2 dEdV1=dEdV1+ eps(1)*ddY1dV1 + eps(2)*dY1dV1 dEdV2=dEdV2+ eps(2)*ddY2dV2 - eps(1)*dY2dV2 где eps(1) и eps(2) находятся функцией
function [eps] = Fr1 = Y(2); eps(1) = dY(1) Fr2 = 1 - Y(1); eps(2) = dY(2) end
В главный модуль генерируемых весов ] Нгуен-Видроу [3]: Hn = power(m,1/n) beta = 0.7*Hn; for i=1:m W1(i)
GetEps(Y,dY)
- Fr1;
- Fr2;
встроена коррекция случайно смещений согласно правилу
beta*W1(i)/abs(W1(i)); beta*W2(i)/abs(W2(i)); -Hn + 2*Hn*rand; -Hn + 2*Hn*rand;
W2(i) B1(i) B2(i) end;
а коэффициент обучения адаптируется в зависимости от текущего значения ошибки: if err>err_old
learnRate = learnRate*0.7; else
learnRate = learnRate*1.05; end;
Эти вставки позволили добиться стабильно успешного выполнения задачи (13) в среднем в пределах 0,5-1 секунды. Графическое сравнение теоретического и нейросетевого решений (рис. 3) доказывает правильность получаемых расчетных данных.
Рисунок 3 - Теоретические (-) и нейросетевые (о) результаты решения задачи (13)
Полученные данные свидетельствуют о высокой производительности разработанного методического комплекса.
ЛИТЕРАТУРА
1. Балабанов А. А. КОПРАС - Техника компьютерного моделирования автоматических систем. Учебный курс. Книга 1. Кишинэу, ТУМ, 2004, 310 с.
2. A. Malek, R. Shekari Beidokhti. Numerical solution for high order differential equations using a hybrid neural network-Optimization method. Applied Mathematics and Computation, 183, 2006, 260-271.
3. Осовский С. Нейронные сети для обработки информации: Пер. с польского. - М.: Финансы и статистика, 2002, 344 с.
4. M. Otadi1, M. Mosleh. Numerical solution of quadratic Riccati differential equation by neural network. Mathematical Sciences Vol. 5, No. 3, 2011, 249-257.
5. M. Shafiee, S. Amani. Optimal control for a class of singular systems using neural network. Iranian Journal of Science & Technology, Transaction B, Engineering, Vol. 29, No. B1, 2005, 3348.
6. Шибанов С.В. Обзор современных методов интеграции данных в информационных системах / Шибанов С.В., Яровая М.В., Шашков Б.Д., Кочегаров И.И., Трусов В.А., Гришко А.К. // Труды международного симпозиума Надежность и качество. 2010. Т. I. С. 292-295.
7. Ерофеев С. А. Синтез оптимального управления сложными конечно-элементными объектами на базе метода Рунге-Кутта в среде MATLAB.
что
д
УДК 621.396.6, 621.8.019.8
Егоров А.М., Новиков П.Г., Кулыгин В.Н.
Национальный исследовательский университет «Высшая школа экономики», Москва, Россия РАЗРАБОТКА ИНТЕРАКТИВНОГО ИНТЕРФЕЙСА ДЛЯ СИСТЕМЫ РАСЧЕТА НАДЕЖНОСТИ И ДОЛГОВЕЧНОСТИ ЭЛЕКТРОННЫХ МОДУЛЕЙ АСОНИКА-К-СЧ
Введение
Система АСОНИКА-К-СЧ представляет собой визуальную среду обеспечения надежности электронных модулей, предназначенную для автоматизации выполнения мероприятий и управления надёжностью радиоэлектронной аппаратуры на этапах ее проектирования, изготовления, эксплуатации и утили-
зации [1]. Система выполнена в виде клиент-серверного приложения, что обеспечивает воспроизводимость расчетов и простоту в эксплуатации [2]. Не смотря на то, что программных средств, аналогичных системе АСОНИКА-К-СЧ, до сих пор не создано ни в России, ни за рубежом, постоянное обновление операционных систем привело к тому,
