Лазарев Ю.Ф. Mатематическое моделирование физических процессов и технических систем в MATLAB

Подождите немного. Документ загружается.

341

Процесс образования блока S-функции состоит из следующих этапов:

1) написания текста S-функции, например, в виде М-файла; текст составляется на основе файла-шаблона

SfunTMPL.m

с учетом заданных уравнений поведения блока;

2) перетаскивания из библиотеки SIMULINK

(

Simulink

►

User-Defined Functions

) стандартного

блока

S-function (рис. 8. 12) в окно блок-схемы, внутри которой будет создаваться новый S-блок;

Рис. 8. 12. Внешний вид заготовки S-функции

3) двойного щелчка мышью на изображении этого блока, что приводит к возникновению на экране

нового окна (рис. 8. 13), содержащего два окошка ввода:

S-function name (имя S-функции), в которое вводится имя новой написанной S-функции;

S-function parameters

(параметры S-функции), в которое вводятся имена или значения тех параметров

блока, которые указаны

в разделе <Параметры> составленного М-файла S-

функции;

Рис. 8. 13. Окно настраивания блока S-function

4) введения в указанные окна имени написанного М-файла S-функции и списка значений параметров S-

функции; если, например, ввести в качестве имени

S_KA

, а в окошко параметров, - строку

J, Ug0, UgSk0

то это окно изменит свой вид на следующий (рис. 8. 14):

Рис. 8. 14. Окно настраивания блока S-function после ввода параметров

342

5) теперь нужно нажать мышью на кнопку <OK> этого окна; если система обнаружит М-файл с

введенным именем в достижимых ей путях, окно (рис. 8. 14) исчезнет, а изображение блока в окне блок-

схемы изменит свой вид (рис. 8. 15) – внутри него возникнет введенное имя S-функции.

S-блок на основе S-функции, которая содержится в М-файле по

имени S_KA.m, создан.

Теперь его можно использовать как обычный S-блок в блок схеме S-модели. При этом на вход этого блока

должен поступать векторный сигнал u. Выходом блока станет векторный сигнал y, который сформирован S-

функцией во внутренней процедуре

mdlOutputs

8.1.6. Примеры создания S-функций

Образуем S-функцию, которая реализует динамические свойства твердого тела в его вращательном движении.

Для описания динамики тела воспользуемся динамическими уравнениями Эйлера в матричной форме:

MωJω

J =⋅⋅×+⋅ )()(

, (8.1)

здесь

- матрица моментов инерции тела относительно осей, связанных с телом; - матрица-столбец из

проекций абсолютной угловой скорости тела на одни и те же оси;

- кососимметричная матрица вида

J ω

)( ×ω

⎥

⎦

⎤

⎢

⎣

⎡

ωω−

ω−ω

ωω−

=×

)(

, (8.2)

составленная из тех же проекций;

- матрица-столбец из проекций вектора момента внешних сил на

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

Поэтому далее иногда тело будет отождествляться с космическим аппаратом.

Создадим М-файл соответствующей S-функции. Назовем его S_DUE.m :

Файл S_DUE.m

function [sys,x0,str,ts] = S_DUE(t,x,M,flag,J,UgSk0)

% S-функция S_DUE Динамических Уравнений Ейлера

% Реализует динамику вращательного движения твердого тела,

% отыскивая вектор абсолютной угловой скорости тела

% по заданному вектору моментов внешних сил,

% действующих на тело

% ВХОД блока:

% M - вектор проекций момента внешних сил на оси

% X, Y i Z связанной с телом системы координат

% ВЫХОД блока:

% y - вектор из шести элементов:

первые три - проекции

% абсолютной угловой скорости om тела на указанные оси,

% последние три - проекции на те же оси

% углового ускорения тела

% Входные ПАРАМЕТРЫ S-функции:

% J - матрица моментов инерции тела в указанных осях;

% UgSk0 - вектор начальных значений проекций

% угловой скорости тела

% Лазарев Ю.Ф., Украина, 18-12-2001

IJ=inv(J); % вычисление обратной матрицы моментов инерции

switch flag,

case 0

[sys,x0,str,ts] = mdlInitializeSizes(UgSk0);

case 1,

sys = mdlDerivatives(t,x,M,J,IJ);

case 3,

sys = mdlOutputs(x,M,J,IJ);

case 9

sys = [];

end

% Конец процедуры

%=======================================================

343

function [sys,x0,str,ts] = mdlInitializeSizes(UgSk0,Ug0)

sizes = simsizes;

sizes.NumContStates = 3;

sizes.NumDiscStates = 0;

sizes.NumOutputs = 6;

sizes.NumInputs = 3;

sizes.DirFeedthrough = 1;

sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

x0 = UgSk0;

str = [];

ts = [0 0];

% Конец процедуры mdlInitializeSizes

%==========================================

function z = mdlDerivatives(t,x,M,J,IJ)

% ВХОДНОЙ вектор M является вектором проекций моментов внешних сил,

% действующих на космический аппарат соответственно по осям X Y Z

% х(1)=om(1); х(2)=om(2); х(3)=om(3);

% z(1)=d(om(1))/dt; z(2)=d(om(2))/dt; z(3)=d(om(3))/dt;

% |Jx Jxy Jxz|

% J = |Jxy Jy Jyz| - матрица моментов инерции КА

% |Jxz Jyz Jz |

om=x(1:3);

z=IJ*(M-cross(om,J*om)); % ДИНАМИЧЕСКИЕ уравнения Ейлера

% Конец процедуры mdlDerivatives

%=================================

function y = mdlOutputs(x,M,J,IJ)

y(1:3)=x;

om=x(1:3);

zom=IJ*(M-cross(om,J*om)) ; % Определения УСКОРЕНИЙ

y(4:6)=zom;

% Конец процедуры mdlOutputs

В качестве входного вектора создаваемого S-блока принят вектор М, состоящий из трех значений текущих

проекций момента внешних сил, которые действуют на тело, на оси системы декартових координат, связанной

с телом. Образуем выходной вектор y из шести элементов: первые три – текущие значения проекций

абсолютной угловой скорости тела, вторые три, - проекции на те

же оси абсолютного углового ускорения тела:

y=[omx,omy,omz,epsx,epsy,epsz].

Создаваемый S-блок рассматривается как непрерывная система (с тремя непрерывными состояниями

x=[omx,omy,omz]

). Поэтому в тексте S-функции изъята процедура mdlUpdate и оставлена процедура

mdlDerivative

, которая фактически является подпрограммой правых частей динамических уравнений Ейлера.

Создадим новое (пустое) окно блок-схемы. Перетянем в него стандартный блок S-функции из раздела

User-

Defined Functions

Дважды щелкнув мышью на изображении этого блока, вызовем его окно настраивания и запишем в него

название М-файла созданной S-функции и его параметры (рис. 7.101). Нажмем кнопку <OK>.

Рис. 8. 15. Окно настраивания блока S_DUE

344

Вследствие этого на изображении блока возникнет надпись имени S-функции, а окно настраивания исчезнет.

Теперь в том же окне с S-блоком S-функции создадим блок-схему для проверки правильности его работы.

Но для этого предварительно нужно продумать условия тестового примера.

Рассмотрим, например, такой случай:

- на тело не действуют внешние силы, т. е

. тело свободно вращается в пространстве;

- оси декартовой системы координат, которая связана жестко с телом, направлены по главным осям

инерции тела; при таком условии матрица моментов инерции будет диагональной;

- тело является динамично симметричным, а его ось фигуры направлена вдоль второй оси (Y) связанной

системы координат; это означает, что матрица моментов

инерции будет иметь вид:

⎥

⎦

⎤

⎢

⎣

⎡

где обозначено:

- экваториальный момент инерции, - момент инерции тела относительно его

оси фигуры (осевой момент инерции тела), причем

в случае, если тело является

сплющенным в направлении оси

;

JJ >

- тело предварительно приведено во вращение с угловой скоростью

вокруг оси его фигуры и имеет

незначительную (по сравнению с

) начальную угловую скорость вокруг оси Х.

Ω

Рассмотрим теоретически движение тела при этих условиях. Уравнения Эйлера в этом случае приобретут

следующий вид:

⎪

⎩

⎪

⎨

⎧

⋅⋅−−=⋅

=⋅

⋅⋅−=⋅

ωω

)(

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

)cos(

Ω⋅=

ωω

; ; , (8.3)

Ω=

)sin(

Ω⋅−=

ωω

причем

−

. (8.4)

Итак, если для образованной модели обеспечить указанные условия, то, если она правильна (адекватна), при

моделировании должны получить результаты соответствующие (8.3).

Добавим в блок-схему блок констант, который формирует нулевой вектор моментов внешних сил (рис. 8. 16), а

также блоки Scope по которым можно проконтролировать результаты моделирования в виде зависимостей от

времени проекций угловой

скорости и углового ускорения тела.

Перед началом моделирования нужно присвоить значение матрицы моментов инерции. Это можно сделать в

командном окне MatLab, вводя строку

J=diag([400,600,400])

Результатом будет возникновение в том же окне записи

J =

400 0 0

0 600 0

0 0 400

345

Рис. 8. 16. Блок-схема проверки работы блока S_DUE

Аналогично нужно ввести вектор начальных условий

UgSk0=[0.001 0.01 0]

Получим

UgSk0 = 1.0000e-003 1.0000e-002 0

Теперь следует перейти к окну блок-схемы, установить параметры интегрирования, указанные на рис. 8. 17 и

запустить блок схему на моделирование. По окончании процесса моделирования, обращаясь к окнам блоков

Scope

, можно убедиться, что созданная модель работает полностью адекватно.

Рис. 8. 17. Параметры процесса интегрирования в S-модели Prov_DUE

Но убедить в этом читателя довольно трудно из-за следующих обстоятельств. Графические результаты блоков

Scope

представлены на черном поле. Поэтому скопировать соответствующее графическое окно означает

получить некачественное изображение на бумаге, вдобавок израсходовав нерациональное количество

красящего материала. Можно скопировать изображение этого окна при помощи команды печати в графическом

окне блока

Scope, но тогда соответствующее изображение займет целый лист, его невозможно уменьшить

средствами текстового редактора. Вдобавок, линии на графике после печати на черно-белом принтере, не будут

отличаться один от другого. На них невозможно нанести надписи, чтобы указать особенности кривых, нельзя

довольно просто изменить стиль линии.

Все это говорит о том, что

наиболее рационально передать результаты в рабочее пространство с помощью

введения выходных портов (рис. 8. 16) и посылания на них тех сигналов, которые нужно представить

графически, а потом построить необходимые графики, пользуясь графическими возможностями MatLab.

Последнее можно сделать непосредственно, пользуясь командами MatLab в его командном окне, но более

рационально сделать это программно, причем желательно объединить в

этой программе все действия:

346

введение значений параметров, начальных условий и т.п.;

- установление параметров интегрирования;

- обращение к S-модели и запуск ее на моделирование;

- обработку полученных результатов моделирования, построение и оформление графиков.

Пример такой программы приведен ниже.

% Prov_DUEupr.m

% Управляющая программа

% для запуска модели Prov_DUE.mdl

% Лазарев Ю.Ф. 18-12-2001

J=[400 0 0; 0 600 0; 0 0 400]; % Ввод значений матрицы инерции

UgSk0=[0.001 0.01 0]; % Ввод начальных значений

% проекций угловой скорости тела

% Установка параметров моделирования

options=simset('Solver','ode45','RelTol',1e-6);

sim('Prov_DUE',5000,options); % МОДЕЛИРОВАНИЕ на S-модели

% Формирование данных и вывод ГРАФИКОВ

tt=tout;

omx=yout(:,1); omy=yout(:,2); omz=yout(:,3);

epsx=yout(:,4); epsy=yout(:,5); epsz=yout(:,6);

subplot(2,1,1)

h=plot(tt,omx,tt,omy,'.',tt,omz,'--');grid

set(h,'LineWidth',2);

set(gca,'LineWidth',2)

set(gca,'FontAngle','italic','Fontsize',16)

title(' Проекции угловых скоростей')

ylabel('радианы в секунду')

legend('omx','omy','omz',0)

subplot(2,1,2)

h=plot(tt,epsx,tt,epsy,'.',tt,epsz,'--');grid

set(h,'LineWidth',2);

set(gca,'LineWidth',2)

set(gca,'FontAngle','italic','Fontsize',16)

title('Проекции угловых ускорений ')

ylabel('1/c2)')

xlabel('Время (с)')

set(gcf,'color','white')

legend('epsx','epsy','epsz',0)

Обратившись к этой программе, получим графики, приведенные на рис. 8. 18. Теперь читатель может наглядно

убедиться в правильности модели.

Рис. 8. 18. Свободное движение симметричного гироскопа (космического аппарата)

Отметим, что промоделированное движение соответствует свободному движению симметричного гироскопа –

его нутационным колебаниям. В гироскоп тело превращает приведение его в сравнительно быстрое вращение

347

вокруг одной из его осей ( ). Матрица моментов инерции принята диагональной, т. е. тело

предполагается динамически сбалансированным относительно осей XYZ. Наконец, моменты инерции

относительно осей, ортогональных оси собственного вращения тела приняты равными. Это означает, что

гироскоп является динамически симметричным телом с осью фигуры Y.

01,0

−

=ω c

8.1.7. Запуск M-программ из S-модели

Следует указать еще один, более удобный, способ объединения S-модели с программами на языке MatLab,

который состоит в возможности вызова М-файлов непосредственно из S-модели специально

предусмотренными для этого средствами.

Например, вызов перед началом загрузки S-модели по имени, к примеру,

MODEL.mdl некоторого М-файла,

например, по имени

PERVdan, который содержит операции присваивания исходных значений всех данных,

можно осуществить, если при создании S-модели с этим именем в командном окне MatLab ввести команду

set_param('MODEL','PreLoadFcn','PERVdan')

Эта команда свяжет файл PERVdan.m c S-моделью MODEL.mdl так, что он будет автоматически

вызываться и исполняться при вызове этой S-модели. Если после выполнения этой команды записать на диск

эту S-модель, то при дальнейших ее вызовах сначала автоматически будет вызван и выполнен файл

PERVdan.m

и лишь после этого на экране возникнет блок-схема S-модели, готовая к моделированию.

Проверить, какой именно m-файл используется в данной S-модели как предварительно выполняемый, можно

путем привлечения команды

get_param('имя S-модели','PreLoadFcn')

Вообще при помощи функции set_param можно установить в S-модели значения многих ее параметров, в том

числе и параметров отдельных ее блоков.

В общем виде обращение к ней может иметь вид

set_param('Имя_S-mодели/Имя_блока','Параметр1',Значение1,

'Параметр2',Значение2,...)

Если указано Имя_блока, то последующие значения присваиваются параметрам этого блока.

Приведем примеры.

Вызов вида

set_param('MODEL','Solver','ode15s','StopTime','3000')

приведет к установлению в S-модели MODEL решателя

ode15s

и времени окончания процесса

моделирования 5000.

Если обратиться к этой функции так:

set_param('MODEL/Уравнение','Gain','1000')

то в блоке Уравнение S-модели MODEL параметру Gain будет присвоено значение 1000.

Команда

set_param('MODEL/Fcn','Position',[50 100 110 120])

установит изображение блока Fcn в S-модели MODEL в прямоугольник в окне бло-схемы с координатами [50

100 110 120].

Если же осуществить такое обращение:

set_param('mymodel/Compute','OpenFcn','my_open_fcn')

то блок

Compute

S-модели

mymodel

будет связан с М-программой MatLab, записанной в файле

my_open_fcn.m

. После этого файл

my_open_fcn.m

будет вызываться на исполнение всякий раз, когда на

изображении блока

Compute дважды щелкнуть мышью.

В случае, когда нужно вызвать некоторый m-файл перед или после проведения собственно моделирования на S-

модели (например, нужно вызвать программу, которая позволяет изменить установленные значения параметров

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

форме), можно установить на свободном месте блок-схемы пустые блоки

Subsystem (из раздела

Ports &

Subsystems

), каждый из которых будет осуществлять вызов соответствующего М-файла.

Эти пустые блоки

Subsystem в поле блок-схемы модели можно соединить с определенной М-программой,

набрав в командном окне команду, аналогичную последней.

348

Попробуем организовать такую форму управления процессом моделирования для предыдущей модели

Prov_DUE

Для этого создадим на основе ранее созданного файла

Prov_DUEupr.m

три отдельных файла

Prov_DUE_Pred. m

осуществляющий предварительное присвоение значений

исходным величинам;

MENU_DUE.m

осуществляющий диалоговое изменение исходных данных;

Graf_DUE.m

осуществляющий вывод результатов моделирования в

графическое окно

Тексты этих программ приведены ниже.

% Prov_DUE_Pred.m

% Программа установки начальных (вшитых) значений

% параметров модели Prov_DUE.mdl

% Лазарев Ю.Ф. 28-01-2004

J=[400 0 0; 0 600 0; 0 0 400]; % Ввод значений матрицы инерции

UgSk0=[0.001 0.01 0]; % Ввод начальных значений

% проекций угловой скорости тела

% Установка параметров моделирования

hi=10; TK=5000;

% MENU_DUE.m

% Программа изменения исходных значений

% параметров модели Prov_DUE.mdl

% Лазарев Ю.Ф. 28-01-2004

k=1;

while k<12

k=menu('Данные для модели Prov_DUE. Что изменить?',...

sprintf('Jx = %g',J(1,1)),...

sprintf('Jy = %g',J(2,2)), ...

sprintf('Jz = %g',J(3,3)),...

sprintf('Jxy = %g',J(1,2)),...

sprintf('Jxz = %g',J(1,3)),...

sprintf('Jyz = %g',J(2,3)),...

sprintf('OMx(0) = %g',UgSk0(1)),...

sprintf('OMy(0) = %g',UgSk0(2)),...

sprintf('OMz(0) = %g',UgSk0(3)),...

sprintf('hi = %g',hi),sprintf('TK = %g',TK),...

'Ничего не менять');

if k==1

J(1,1)=input([sprintf('Текущее значение J(1,1)=%g;

',J(1,1)),'Установите новое значение J(1,1)=']);

end

if k==2

J(2,2)=input([sprintf('Текущее значение J(2,2)=%g;

',J(2,2)),'Установите новое значение J(2,2)=']);

end

if k==3

J(3,3)=input([sprintf('Текущее значение J(3,3)=%g;

',J(3,3)),'Установите новое значение J(3,3)=']);

end

if k==4

J(1,2)=input([sprintf('Текущее значение J(1,2)=%g;

',J(1,2)),'Установите новое значение J(1,2)=']);

end

if k==5

J(1,3)=input([sprintf('Текущее значение J(1,3)=%g;

',J(1,3)),'Установите новое значение J(1,3)=']);

end

if k==6

J(2,3)=input([sprintf('Текущее значение J(2,3)=%g;

',J(2,3)),'Установите новое значение J(2,3)=']);

end

if k==7

349

UgSk0(1)=input([sprintf('Текущее значение OMx(0)=%g;

',UgSk0(1)),'Установите новое значение OMx(0)=']);

end

if k==8

UgSk0(2)=input([sprintf('Текущее значение OMy(0)=%g;

',UgSk0(2)),'Установите новое значение OMy(0)=']);

end

if k==9

UgSk0(3)=input([sprintf('Текущее значение OMz(0)=%g;

',UgSk0(3)),'Установите новое значение OMz(0)=']);

end

if k==10

hi=input([sprintf('Текущее значение hi=%g; ',hi),...

'Установите новое значение hi=']);

end

if k==11

TK=input([sprintf('Текущее значение TK=%g; ',TK),...

'Установите новое значение TK=']);

end

J(2,1)=J(1,2); J(3,1)=J(1,3); J(3,2)=J(2,3);

% Graf_DUE.m

% Программа построения в графическом окне графиков

% результатов работы модели Prov_DUE.mdl

% Лазарев Ю.Ф. 28-01-2004

% Формирование данных и вывод ГРАФИКОВ

tt=tout;

omx=yout(:,1); omy=yout(:,2); omz=yout(:,3);

epsx=yout(:,4); epsy=yout(:,5); epsz=yout(:,6);

StrJ=[sprintf('Jx= %g; ',J(1,1)),sprintf('Jy= %g; ',J(2,2)),sprintf('Jz= %g; ',J(3,3))]

StrU=[sprintf('OMx= %g; ',UgSk0(1)),sprintf('OMy= %g; ',UgSk0(2)),sprintf('OMz= %g;

',UgSk0(3))]

StrJ1=[sprintf('Jxy= %g; ',J(1,2)),sprintf('Jxz= %g; ',J(1,3)),sprintf('Jyz= %g;

',J(2,3))]

subplot(2,1,1)

h=plot(tt,omx,tt,omy,'+',tt,omz,'--');grid

set(h,'LineWidth',2);

set(gca,'LineWidth',2)

set(gca,'FontAngle','italic','Fontsize',16)

title(' Проекции угловых скоростей')

set(gca,'Fontsize',14)

ylabel('радианы в секунду')

legend('omx','omy','omz',0)

xlabel([StrJ,' ',StrJ1,' ',StrU])

subplot(2,1,2)

h=plot(tt,epsx,tt,epsy,'+',tt,epsz,'--');grid

set(h,'LineWidth',2);

set(gca,'LineWidth',2)

set(gca,'FontAngle','italic','Fontsize',16)

title('Проекции угловых ускорений ')

set(gca,'Fontsize',14)

ylabel('радианы/c^2')

xlabel(' Время (с) ')

set(gcf,'color','white')

legend('epsx','epsy','epsz',0)

Теперь введем в блок-схему Prov_DUE два блока Subsystem (рис. 8. 19).

350

Рис. 8. 19. Блок-схема S-модели Prov_DUE1

Произведем в этих блоках следующие изменения:

1) изменим названия блоков (подписи под их изображениями); под первым блоком запишем MENU, под

вторым GRAFIKI;

2) дважды щелкнув на каждом из блоков, откроем их блок-схемы и уничтожим все их содержимое, -

сделаем блоки пустыми;

3) выйдя вновь в окно блок-схемы, щелкнем на

изображении блоков правой клавишей мыши; при этом

появится список команд, из которых следует выбрать команду Mask Subsystem; при этом появится окно

Mask Editor (рис. 8. 20);

4) введем в поле Drawing commands вкладки Icon этого окна команду

disp с указанием в качестве ее

аргумента текста, который мы желаем поместить на изображении блока; для первого блока запишем

Исходные данные, для второго – Вывести графики;

5) закроем окно

Mask Editor

, нажав в нем кнопку

OK.

В результате получим изображения блоков, представленные на рис. 8. 19.

Рис. 8. 20. Окно Mask Editor

Введем в командном окне команды, связывающие составленные М-программы с S-моделью: