Моделируем распространение коронавирусной инфекции COVID-19 в Simulink - часть 2

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

В этой статье мы разберем, как устроена модель (спойлер: очень интересно с инженерной точки зрения) и освоим на ее примере несколько мощных нововведений последних релизов MATLAB.

Изначально эту модель я увидел в блоге разработчика из Матворкс Гая Ролё (Guy Rouleau). Меня заинтриговало применение конечных автоматов Stateflow для решения это задачи, в частности ряда хитростей, которые сильно упрощают построение модели, но до которых нужно еще додуматься.

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

Напоминаю, что вы можете самостоятельно поиграться с этой моделью, для чего можете оперативно получить от нас пробную версию MATLAB R2020a на время карантина.

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

Модель движения популяции

По сути, все симуляция происходит внутри Stateflow-диаграммы Propagation Simulator, которая представлена на картинке. При этом диаграмма имеет всего одно состояние Integ, в котором симулятор находится все время.

Обратите внимание, что Integ - это не обычное состояние Stateflow, это Simulink-состояние, которые появились в Stateflow несколько релизов назад. Внутри этого состояния находится Simulink-диаграмма, которая пересчитывается на каждом такте работы модели (поскольку автомат его не покидает).

Внутри состояния Integ находится нехитрая схема из двух интеграторов, которая моделирует дифференциальное уравнение движения твердого тела.

Обратите внимание, что у каждого интегратора сверху есть странная иконка с подписью. Это означает, что мы имеем доступ к состояниям этих интеграторов (к их накопленным значениям). Соостояние первого интегратора обозначено в его настройках как v0 (скорости [Vx,Vy]), второго - pos0 (положение [x,y]). Этим доступом к состояниям мы воспользуемся дальше (это самая большая хитрость).

Поскольку наши люди-шарики перемещаются на плоскости, то каждый интегратор просчитывает одновременно и координату x и координату y. Таким образом это уравннеие Simulink векторизовано, посколько в каждый момент времени просчитывает вектор [x,y], что позволяет нам решать только одно дифференциальное уравнение, а не два одинаковых. Получается, что мы используем матричные возможности MATLAB/Simulink, чтобы упростить модель и повысить ее быстродействие.

Давайте еще раз посмотрим на анимацию движения людей, которую разбирали в прошлый раз.
(Нажмите на картинку, чтобы анимировать).

Заметьте, что точки не просто движутся на плоскости, они еще и сталкиваются и отскакивают друг от друга и от стен. Чтобы такое просимулировать, нам нужно на каждом такте работы модели не только просчитывать координаты всех точек (чем занимается состояние Integ), но и определять столкновения и просчитывать изменение скорости каждой точки при контакте.

За это отвечают специальные функции, созданные в этой же диаграмме Statelow, вы можете увидеть их в левой части диаграммы:

• detectHit - MATLAB-функция, которая определяет, что призошло столкновение точек друг с другом
• contact - Simulink-функция, которая рассчитывает изменение движения точек при контакте
• detectHitWall - MATLAB-функция, которая определяет, что призошло столкновение точек со стеной
• contactWall - MATLAB-функция, которая рассчитывает движение точки при контакте со стеной

Вызываются все эти функции на каждом такте работы модели, потому что их вызов прописан на переходах состояния Integ из себя в себя же. Это такая небольшая хитрость, чтобы заставить эти функции вычисляться на каждом такте.

При этом в каждую функцию передаются матрицы pos и v, у каждой из которых в нашем случае размерность [100,2]. То есть они содержат информацию о движении на плоскости xy для каждой из ста точек, которые мы просчитываем. Снова векторизация в действии! Просчитываем все точки за раз.

Дальше уже сама функция внутри себя решает, как обрабатывать эти матрицы точек. Например, внутри MATLAB-функции detectHit реализованы два вложенных цикла for, которые перебирают все точки и вычисляют, есть ли контакт между ними.

Схожим, но более сложным образом устроена Simulink-функция contact, внутри которой можно найти сразу 3 цикла, реализованных с помощью блоков For Iterator Subsystem. Рекомендую поиследовать эту функцию самостоятельно, чтобы расширить свое владение подобными конструкциями.

Продравшись через все итераторы и циклы в сердце диаграммы contact, мы находим математику, которая рассчитывает конакт твердых тел.

Здесь все по-честному - происходит преобразование систем координат, вычисление нормальных и тангенциальных скоростей, пересчет скоростей при контакте. Если доберетесь до сюда, то можете переключить вариант подсистемы Contact Dynamic с Simple Model на Physical Model, в которой реализован закон сохранение импульса и потери энергии при ударе. В такой конфигурации шарики будут чуть замедляться при каждом столкновении и выбивать друг друга с места. Интересно, как это повлияет на результат распространения?

Кстати, с вариантом Physical Model эту модель можно использовать не только для моделирования распространения инфекции, но и для симуляции бильярда, кёрлинга и даже сыпучих сред! Так что можем с вами придумать еще что-нибудь интересное на основе этой модели, предлагайте свои идеи в комментариях.

Так вот, вернемся к интеграторам Integ. На втором переходе из этого состояния написано действие:

[Integ.v0, b1, b2] = contact(pos,v);

То есть один из результатов функции contact мы записываем напрямую в состояние интегратора Integ.v0. Иначе говоря, при каждом контакте мы пересчитываем скорости всех точек с помощью функции contact, обновляем их внутри нашего дифференциального уравнения в Integ, а координаты точек автоматически высчитываются вторым интегратором pos0.

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

Помимо гибкости есть еще одно достоинство такого подхода - это оптимизация быстродействия модели. Да, при моделировании контактов мы вынуждены использовать вложенные циклы, без них задачу решить намного сложней, ведь нам надо отслеживать контакты всех точек со всеми. Но, что касается моделирования динамики, то она одинаковая для каждой точки, и здесь вместо цикла for на сто итераций мы ипользуем одно-единственное уравнение в Integ, в которое в качестве скорости мы записываем матрицу [100,2] и одновременно просчитываем 200 уравнений, что работает намного быстрее циклов благодаря матричным возможностям MATLAB.

Отслеживание статуса заражения

Мы с вами разобрали интересную реализацию модели движения большого количества объектов. Осталось чуть-чуть - рассмотреть, как моделируются процессы заражения и выздоровления.

Статус заражения всех людей хранится в векторе status, который состоит из 100 элементов и обновляется на каждом такте расчета модели с помощью графической функции updateStatus, в которой реализован поток управления Stateflow.

Тут все просто - если зараженный (status: 1) сталкивается со здоровеньким (status: 0), то происходит заражение последнего и его статус меняется на 1.

А на последнем переходе этого потока управления вызывается Simulink-функция recover, описанная на верхнем уровне диаграммы.

Тут я тоже применил небольшую хитрость. Время выздоровления контролируется с помощью интегратора. Как только человек становится зараженным, на входе интегратора появляется единица, и он начинает копить значение со скоростью 1 ед/сек. Как только накопленное значение становится больше времени выздоровления, встроенная MATLAB-функция меняет статус человека с 1 на 2 (выздоровел и имеет иммунитет). И конечно же функция recover тоже векторизована и обрабатывает статус всех 100 человек одновременно.

Заключение

На этом мы заканчиваем рассмотрение модели Simulink, надеюсь вам было интересно и вы узнали что-то новое. Так что скачивайте модель и экспериментируйте. Если что непонятно, спрашивайте в комментах. В следующей части рассмотрим, как устроена визуализация распространения и анимация графиков.

P.S. Если вы так же любите Stateflow, как и я, смотрите свежий вебинар

Вебинар «Разработка конечных автоматов и управляющей логики и их запуск на микроконтроллерах и ПЛИС»