4.7. Ускорение вычислений
Этот раздел описывает методы повышения быстродействия при выполнении программы и управление памятью:
векторизация циклов;
предварительное размещение векторов.
Под векторизацией понимается преобразование циклов for и while к эквивалентным векторным или матричным выражениям. При векторизации алгоритма ускоряется выполнение m-файла.
Пример 4.4. Табулирование функции. Одним из способов вычисления 1001 значение функции синуса на интервале [0 10] может быть реализован с помощью оператора цикла:
i = 0;
for t = 0:.01:10
i = i + 1;
y(i) = sin(t);
end
Второй способ: применение векторизованной формы
t = 0:.01:10;
y = sin(t);
Именно второй способ вычисления выполняются намного быстрее, и такой подход в системе Matlab является предпочтительным.
В системе Matlab есть возможность для существенного сокращения времени выполнения программы за счёт предварительного размещения массивов для выходных данных. Для этого необходимо использовать операции zeros, struct и cell. Предварительное выделение памяти позволяет избежать фрагментации памяти при работе с большими матрицами. Предварительное выделение памяти позволяет определить непрерывную область, достаточную для проведения всех вычислений.
Существует несколько подходов к повышению эффективности использования памяти, они рассмотрены ниже. В системе Matlab предусмотрено пять функций для работы с памяти:
команда clear - удаление переменных из оперативной памяти;
команда pack - запись текущих переменных на диск и последующей их загрузкой;
функция quit - по мере необходимости выход из системы Matlab с освобождением всей памяти;
команда save - сохранение переменных в файле;
команда load - считывание данных из файла.
Команды save и load работают быстрее, чем команды ввода-вывода.
4.8. Примеры программирования с использованием средств пакета Matlab
Пример 4.5. Моделирование дискретного фильтра Калмана. Фильтр Калмана используется для построения оценок процессов при косвенных измерениях с ошибками. Подробно алгоритм описан в [14].
Текст m-сценария имеет вид:
%Дискретный фильтр Калмана
A=[0 1;0.05 0.9]; % матрица динамики
S=[0 1]; % матрица канала измерений
Q=[0.1 0;0 0.325]; V=0.09; % интенсивности
P=[0.03 0;0 0.05]; %начальная дисперсия
E=eye(2); %единичная матрица
x=[5 7.5]'; % начальные условия для объекта
xx=[0 0]'; % начальные условия для оценок
T=[]; X=[]; Y=[]; %объявление переменных
i=0;
while i<=20
T=[T i]; X=[X x]; Y=[Y xx];%формирование результата
x=A*x+Q*randn(2,1); % модель объекта
y=S*x+V*randn; %модель канала наблюдений
P1=A*P*A'+Q*Q'; %дисперсия прогноза
K=P1*S'*inv(S*P1*S'+V*V); %коэффициенты передачи фильтра
P=(E-K*S)*P1;% дисперсия оценки на очередном шаге
xx=A*xx+K*(y-S*A*xx);%вычисление оценки
i=i+1;
end;
subplot(2,1,1);
plot(T,X(1,:)',T,Y(1,:)');%графики 1-ой компоненты и ее оценки
grid;
subplot(2,1,2);
plot(T,X(2,:)',T,Y(2,:)');%графики 2-ой компоненты и ее оценки
grid;
Результаты вычислений приведены на рис 4.2; из них видно, что оценки процессов начинают удерживаться вблизи истинных значений практически, начиная с 1-го и 2-го такта.
Рис. 4.1 - Результаты моделирования
Пример 4.6. Модальное управление. Смысл модального управления заключается в следующем. Необходимо для линейной динамической системы
,
(4.1)
найти
управление в виде
,
такое, чтобы матрица замкнутой динамики
имела заданные желаемые собственные
числа. Эта задача решается средствами
пакета Matlab.
Выбор матрицы
осуществляется с помощью m-функции
пакета с именем place.
Кроме того, программа вычисляет матрицу
управляемости Калмана с помощью функции
ctrb (m-файл,
реализующий эту функцию приведен в
примере 4.1).
Ранг матрицы
управляемости вычисляется с помощью
функции rank
(исходная система является управляемой,
если ранг равен размерности исходного
объекта). Решение дифференциального
уравнения замкнутой системы
(4.2)
осуществляется с помощью функции ode45, реализующей метод Рунге-Кутта численного интегрирования дифференциальных уравнений (более подробно проблема численного решения дифференциальных уравнений описана в [3]). Ниже приведен текст m-сценария:
%Модальное управление
A=[0 1; 2 3];
B=[4 3]';
c=ctrb(A,B); %вычисление матрицы управляемости
rnk=rank(c);%вычисление ее ранга
disp('Ранг=')
disp(rnk);
disp('Если ранг равен 2 то система управляема')
r=[-0.3+i -0.3-i]';%желаемые собственные числа
K=place(A,B,r); %вычисление матрицы К
Phi=A-B*K;
rr=eig(Phi); %вычисление собственных чисел матрицы Phi
disp('Собств. числа=')
disp(rr) % вектор rr должен совпасть с r
x0=[1 0.1];
interval=[0 10];%временной интервал моделирования
[T,X]=ode45('model',interval,x0,[],Phi);
subplot(1,2,1)
plot(T,X)%графики переходных процессов
grid
subplot(1,2,2)
plot(X(:,1),X(:,2)) %фазовый портрет
grid
Используемая программой ode45 m-функция model, реализующая правую часть дифференциального уравнения (4.2), имеет вид:
function f=model(t,x,flag,L);
f=L*x;
Здесь t - независимая переменная; x - зависимая переменная; flag - строковая переменная, которая имеет вид возвращаемой информации; L - дополнительный параметр (их может быть несколько, разделяются запятыми), в m-сценарии это матрица Phi.
На рис 4.3 приведены переходные процессы в системе и фазовый портрет. Так как желаемые собственные числа имели отрицательные действительные части, то видим, что замкнутая система устойчива.
Рис. 4.2 - Переходные процессы и фазовый портрет
Пример
4.7.
Оптимальное управление по квадратичному
интегральному критерию. Предполагается,
что модель объекта описывается уравнением
(4.1). Требуется найти управление
,
минимизирующее квадратичный интегральный
критерия вида:
,
где
и
- весовые матрицы. Оптимальное управление
имеет вид:
,
где
матрица
вычислятся с помощью функции пакета
Matlab
lqr.
Входные переменные этой функции являются
матрицы
,
,
и
.
Ниже приведен текст m-сценария:
%Оптимальное управление
global Phi
A=[0 1; 2 3]; B=[4 3]';
С =[1 0; 0 1]; D=1;
K=lqr(A,B,C,D)% вычисление оптимальной матрицы K
Phi=A-B*K;
eig(Phi) %вычисление собственных чисел
x0=[1 1];%начальные условия
tt=[0 3]; %временной интервал моделирования
[T,X]=ode45('fun',tt,x0);
plot(T,X)%графики переходных процессов
Используемый программой ode45 m-функция fun, реализующая правую часть дифференциального уравнения (4.2) имеет вид:
function f=fun(t,x);
global Phi
f=Phi*x;
return
Отметим, что в этом примере матрица Phi объявлена глобальной переменной. На рис 4.4 приведены переходные процессы.
Рис. 4.3 - Переходные процессы в оптимальной системе
Пример
4.8. Минимизации
без ограничений. Рассмотрим
задачу нахождения значений переменных
и
,
обеспечивающих решение задачи минимизации
функции:
.
Нахождение решения производится в соответствии со следующими этапами.
1. Составление m-файла (с именем fun), реализующего вычисление значения целевой функции:
function f=fun(x)
f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);
2. Составление программы с использованием подходящей функции пакета (в данном случае - функции fminunc):
x0=[-1,1];% начальные значения
options=optimset('LargeScale','off');
[x,fval,exitflag,output] = fminunc(‘fun’,x0,options)
Выполнение программы после нескольких шагов приведет к следующему результату:
Optimization terminated successfully
Current search direction is a descent direction, and magnitude of directional derivative in search direction less than 2*options TolFun
x =
0.5 -1.0000
Значение функции в точке х возвращается через fval:
fval=1.3030e-10,
exitflag > 0.
Значение exitflag=l (exitflag > 0) дает информацию, что найдена точка минимума (возможно, локального). Выходная структура (информация о результатах оптимизации) определяется идентификатором output:
число выполненных итераций (iterations): 7;
число вычислений функции (funcCount): 40;
шаг поиска (stepsize): 1;
степень оптимальности найденного решения (firstorderopt) - норма вектора-градиента в точке найденного решения:
8.1995е-004.
С помощью выбора переменных в опциях возможно в fminunc изменить характеристики оптимизационного алгоритма, например как x = fminunc(fun,x0,options). Структура options включает в себя значения точности окончания расчета и выбор алгоритма. Структура options может быть задана с помощью функции optimset: options = optimset('LargeScale','off'); В данном примере исключается принимаемый по умолчанию алгоритм для задачи большой размерности и был использован алгоритм средней размерности. Для задач большой размерности используется функция: options = optimset('LargeScale','on')
Пример 4.9. Минимизация с ограничениями. Пусть оптимизируемая функция будет таже, что и в примере 4.9. Ограничения имеют вид:
Определим функцию ограничений confun с ограничениями в форме неравенств:
function [c,ceq]=confun(x)
%огр. в форме неравенств
c=[1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10];
ceq=[];%в виде равенств ограничений нет
Основная программа:
x0=[-1,1];%начальные значения
options=optimset('LargeScale','off');
[x,fval]=fmincon('fun',x0,[],[],[],[],[],[],'confun',options)
Если минимизация осуществляется с дополнительным диапазоном в векторной форме:
,
то основная программа изменится и будет следующей:
x0=[-1,1];%начальные значения
lb=[0,0];%нижние границы
ub=[];%верхних ограничений нет
options=optimset('LargeScale','off');
[x,fval]=fmincon('fun',x0,[],[],[],[],lb,ub,'confun',options)
В последней строке четыре знака [],[],[],[] означают, что в задаче отсутствуют линейные ограничения в векторно-матричном форме (в виде равенств и неравенств) типа (4.3) и (4.5).
Пример 4.10. Метод наименьших квадратов (МНК) для переопределенной системы
с ограничениями в векторно-матричной форме:
(4.3)
(4.4)
Пример программы для этой задачи с исходными данными имеет вид:
C=[0.951 0.762 0.615
0.231 0.45 0.79
0.61 0.185 0.922
0.48 0.82 0.738];
d=[0.578, 0.353,0.8131,0.0098]';
A=[0.203 0.272 0.747
0.198 0.199 0.445];
b=[0.5251,0.2026]';
lb=-0.1*ones(3,1);
ub=2*ones(3,1);
x=lsqlin(C,d,A,b,[],[],lb,ub)
Cимволы [] и [] в последней строке означают, что в задаче отсутствуют линейные ограничения в виде равенств вида
,
(4.5)
где
- матрица,
- вектор.
Пример 4.11. Определение минимакса. Постановка задачи формулируется следующим образом:
.
Программа для этой задачи в случае 5-ти мерной вектор-функции имеет вид:
x0=[0.1; 0.1];% начальные условия
[x, fval]=fminimax('myfun',x0)
Программа m-файла для вектор-функции цели:
function f=myfun(x)
f(1)=2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304;
f(2)=-x(1)^2-3*x(2)^2;
f(3)=x(1)+3*x(2)-18;
f(4)=-x(1)-x(2);
f(5)=x(1)+x(2)-8;
В задачах на минимакс могут учитываться как линейные так и нелинейные ограничения в виде равенств и неравенств.
Пример 4.12. Линейное программирование. Постановка задачи следующая:
(4.6)
при
дополнительных ограничениях вида (4.3)
– (4.5). В (4.6)
- вектор-столбец,
- вектор-столбец.
Конкретная программа при отсутствии линейных ограничений в виде равенств имеет вид:
f=[-5; -4; -6];
A=[1 -1 1;3 2 4;3 2 0];
b=[20; 42; 30];
lb=zeros(3,1);
ub=[];
x=linprog(f,A,b,[],[],lb,ub)
Пример 4.13. Квадратичное программирование. Постановка задачи следующая:
(4.7)
при
дополнительных ограничениях вида (4.3)
– (4.5). В (4.7)
- матрица,
- вектор-столбец,
- вектор-столбец.
Ниже приведен пример программы.
H=[1 -1; -1 2];
f=[-2; -6];
A=[1 1; -1 2; 2 1];
b=[2; 2; 3];
x=quadprog(H,f,A,b)
В этом примере присутствуют только линейные ограничения в форме неравенств.
Отметим, что задачи максимизации функций сводятся к задачам минимизации тех же функций, но с противоположным знаком.
Пример
4.14.
Оптимизация
нормы матрицы при ограничениях в виде
линейного матричного неравенства (типа
неравенства Ляпунова). Необходимо найти
матрицу
как результат решения следующей задачи:
при ограничениях:
,
где
- заданные матрицы. Пример программы
решения подобной задачи имеет вид:
Q=[-1 -2; -2 -1];
P=[1 -1; -1 -3];
M=[10 0; 0 10];
X=basiclmi(M,P,Q,'Xmin')
Пример 4.15. Оптимизация следа матрицы при ограничениях в виде квадратичного матричного неравенства (типа неравенства Риккати). Необходимо найти матрицу как результат решения следующей оптимизационной задачи:
при ограничениях:
,
где
- заданные матрицы. Пример программы
решения данной оптимизационной задачи
имеет вид:
A=[-1 -2 1; 3 2 1;1 -2 -1];
B=[ 1 0 1]’;
Q=[1 -1 0; -1 -3 -12; 0 -12 -36];
setlmis([])
X=lmivar(1,[3 1])
lmiterm([1 1 1 X],1,A,'s')
lmiterm([1 1 1 0],Q)
lmiterm([1 2 2 0],-1)
lmiterm([1 2 1 X],B',1)
LMIs=getlmis
C=mat2dec(LMIs,eye(3))
option=[1e-5,0,0,0,0]
[copt,xopt]=mincx(LMIs,C,option)
Xopt = dec2mat(LMIs,xopt,X)%фомирование результата
%в виде матрицы
Отметим, что пояснения смысла операторов, используемых в данной программе можно найти в комментариях, выполнив в командном окне команду help «имя оператора».
Пример 4.16. Построение аппроксимирующей нейронной сети (обобщенно-регрессионной). Ниже приведен пример программы с конкретными исходными данными..
in=[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15]%табличные значения
%аргумента
d=[0.2 3 2 5 6 7 8 5 4 3 2 1 1 2 2 ]%табличные значения
%функции
nf=newgrnn(in,d,0.4);%настройка нейронной сети
%по заданной обучающей выборке
%0.4 – допустимая ошибка
x=0:0.1:15 %входной массив
y=sim(nf,x);%опрос (моделирование) сети
plot(x,y,in,d,'*') %построение графика
На рис. 4.5 приведен график аппроксимирующей функции, построенной с помощью нейройной сети.
Рис. 4.4 - График аппроксимирующей функции
Пример 4.17. Линейная прогнозирующая нейронная сеть. Пример программы:
T=[14.0 14.7 14.7 15.1 15.2 15.3 15.8 16.0 17.2 ... %значения
20 20.2 23.1 23 22 20.9 22.1 22.2 24.1 28 27.9 ...%временного
26.2 24.5 22 20 22.2 23.3 22.6 21.4 22.8 22.2 ...%ряда
24.3 28.4 28.7 26.6 21.6 21.1 20.5 23.2 24.7 25.5];
Q=length(T);
time=1:Q;
M=max(T);
T=T./M; %масштабирование исходного временного ряда
P=zeros(5,Q);%5 - глубина памяти алгоритма прогноза
P(1,2:Q)=T(1,1:(Q-1));
P(2,3:Q)=T(1,1:(Q-2));
P(3,4:Q)=T(1,1:(Q-3));
P(4,5:Q)=T(1,1:(Q-4));
P(5,6:Q)=T(1,1:(Q-5));
net=newlind(P,T);%создание нейронной сети
%(проектирование нового линейного слоя)
a=sim(net,P);%опрос сети
e=T-a;%вычисление ошибки
subplot(2,1,1)
plot(time,M*a,time,M*T,'*')%построение графика
%прогнозируемых значений
subplot(2,1,2)
plot(time,M*e,time,0)%построение графика ошибок прогноза
На рис. 4.6 в верхнем окне приведен график значений прогноза и символом «*» обозначениы значения временного ряда, в нижнем окне приведен график ошибки прогноза
Рис. 4.5 - График значений прогноза и ошибок прогноза
Подробней вопросы построения нейронных сетей отражены в учебнике [1].