Мой любимый автор книг по математике Очков Валерий Федорович предложил задачу для среды SimInTech. На самом деле он вызвал меня на дуэль и предложил выяснить кто лучше сделает решение определенной задачи в виде структурной схемы. 

Надо сказать, что структурное моделирование предназначено в первую очередь для моделирования процессов во времени, во вторую — для создания алгоритмов управления сложными техническими объектами, а все остальные математические задачки — это к MathStudio или Wolfram. Структурное моделирование в виде диаграмм, разумеется, позволяет решить любую задачу, но, если задача не динамическая, это уже будет извращение.

С другой стороны, у нас есть примеры, когда среда структурного моделирования использовался как замена Coreal Draw в проектах для авиации, а я сам лично писал статью как можно решать с помощью структурного моделирования решать задачу по построению оптимальных маршрутов муравьев, которую обычно решают AnyLogic. Да и существующий в пакете SimInTech модуль расчёта надежности с графической подготовкой данных, то же выглядит не как нативное приложение для среды моделирования в виде структурных схем. 

Поэтому мы извращений не боимся и смело беремся за операцию на глазном дне через задний проход. Стрелять так стрелять!

В этой статье мы решим задачу в системе динамического моделирования двумя способами: 

1 «Классика» – динамической модель в виде структурной схемы.

2 «Изврат» - метод стрельбы по площадям.

Сама задача и решение выглядит так: 

Есть фигура (см. рис. 1), нужно определить центр масс.

Классика. Как это работает?

Расчет площади фигуры — это совсем не к динамическому моделированию. Мы исследуем процессы, в которых на основании данных по скорости мы получаем графики изменения параметров процесса по времени. Для этого мы решаем дифференциальные уравнения. Казалось бы, как связаны площадь фигур и дифференциальные уравнения? А связаны они через интеграл. В самом деле, центр масс связан с площадью, а площадь фигуры под кривой может быть рассчитана как определенный интеграл функции, описывающий кривую (см. рис. 2):

Зная функцию кривой, можно выполнить расчет площади под этой кривой с помощью интегрирования. А в моделировании динамических систем мы только этим и занимаемся: постоянно интегрируем всё, что движется, а что не движется — двигаем и тоже интегрируем. Поскольку в исходной задаче нужно насчитать площадь фигуры на основе круга, то нам достаточно составить уравнение, в котором какая-либо величина будет меняться по круговой траектории.  Как известно, площадь круга отвечает уравнению:

А значит можно записать выражение для Y:

Если принять, что Х = t — это время, тогда можно создать модель, которая нарисует нам часть окружности. В условиях задачи R = 1 м = 1000 мм. Пусть будет 1000 мм, тогда, чтобы нарисовать круг, время моделирования у нас тоже должно быть 1000 сек. Для людей, которые слушали физику в школе и привыкли к физическим размерностям, это немного странно — рисовать круг, когда по одной оси время, по другой — расстояние. Но математика — наука абстрактная, и тут можно всё.

Зададим константу, которая будет у нас радиусом в квадрате: R² = 1000² = 1E06, используя источник времени, возводим это время в квадрат и отнимает от константы, после этого осталось вычислить корень и все – четвертинка круга нарисована.

Теперь, если добавить на схему интегратор, то он в конце расчета на выходе будет выдавать нам площадь этой самой четвертинки, поскольку у нас время 1000 секунд и радиус = 1000, то площадь будет мм². 

Если бы мой учитель физики в советском техникуме это увидел я бы получил подзатыльник, по оси Х - расстояние по оси Y - время, а площадь это мм². Но мы работаем с математикой, где нет размерности, и мы можем легко это понять, простить и принять.

Что бы не мучиться с сотнями тысяч миллиметров, в интеграторе я добавил перевод мм² в м², поделив все на 1e6 в коэффициенте интегратора, у меня получилась такая схема:

Как видим, при расчете методом Эйлера с шагом 0.1 секунды или 0.1 мм мы получили результат с приемлемой погрешностью в четвертом знаке после запятой. Для сравнения я добавил на схему константу, где рассчитал площадь четвертинки круга с радиусом 1, которая равна:  

Зная площадь, можно посчитать координату центра масс по той оси, которая у нас используется для интегрирования, ось х. Формула центра масс, известная из школы:

Интеграл у нас уже есть, поэтому такую формулу можно реализовать с минимальными доработками в схеме:

У меня получилось 0.424. Чтобы проверить, насколько наша теория соответствует суровой правде жизни, посчитаем центр масс, который нам известен. Например, если бы у нас была не четверть окружности, а половина окружности, то центр массы находился бы в точке 0.

Назад в прошлое

Но там, где четверть круга, там и половина. Из графика на рисунке 4 очевидно, что если бы начало графика было не в 0, а в точке -1000 сек, то по уравнению окружности мы бы как раз получили бы полкруга. То есть нам нужно забраться в прошлое на 1000 секунд и начать считать с отрицательного времени. В жизни такого не бывает, только в голливудском кино.

Но у нас SimInTech, и здесь, как в Голливуде, возможно всё, например, можно поменять время начала расчета в настройках и стартовать не с 0, а из прошлого.

В этом случае на графике линии вместо четвертинки у меня получится идеальная полуокружность, а график центра масс приходит в 0, что и требовалось доказать. Расчет центра масс происходит правильно. Аллилуйя!

Магия и ее разоблачение

Для тех, кто не понял, почему изменяется центр масс, я поясню. В отличие от стационарных расчетов, где окружность готовая и неизменная, в нашей системе она строится на во время моделирования процесса, прямо по времени «на ходу». Мы с каждым шагом расчета добавляем точки и каждый раз считаем новое положение центра масс получившейся фигуры и отображаем это положение на графике 9.

В начальный момент времени у нас вместо круга только одна точка, при начальном значении времени ее координата -1000 секунд, что соответствует координате – 1000 мм или – 1 м (мы пересчитали мм в м в интеграторе). Там же находится и центр масс в начале расчета. Постепенно в график добавляются точки, и он движется по окружности, и центр масс сдвигается. В итоге, когда график заканчивается идеальным симметричным полукругом, центр масс перемещается в 0.

Так с помощью магии SimInTech время превращается в расстояние, не говорите только это своему учителю физики. 

Кстати, именно потому, что задача симметричная мы можем не считать полный круг. В самой задачи тоже требуется посчитать положение центра масса по осу Y. см. рисунок 1.

Поскольку в условиях задачи угол 180 градусов — это как раз полкруга, мы его развернём на 90%, и наша ось X, которая время T, превратится в ось Y. Найдя центр масс по этой оси Х в схеме, мы и получим нужную нам величину Y.

Время у нас превратилось в Х, а Х в Y. Кручу, верчу,  запутать хочу!

Начинаем резать

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

Формула ровно такая же, только меняется радиус, в итоге у меня получается вот такая простая схема расчета:

К исходной схема расчет круга R= 1 м добавляется схема расчета круга R=0.9 мм

Потом от первого круга отнимается второй, и всё это отправляется на два интегратора, рассчитывающих центр масс по формуле. В результате расчета на выходе мы получаем занесение 0.60532671, что практически совпадает с результатом, посчитанным по теоретическим формулам. 0.60534722.

А в результате на графике мы получили две дуги, которые и ограничивают нашу фигуру.

Если раньше это решение должно возмущать учителя физики из-за произвола с размерностями величин, то сейчас должен обплеваться учитель математики, поскольку вторая окружность явно нарушает законы математики. Обратите внимание на график после 900 секунды, например 901 секунда, в этом случае, согласно формуле, мы должны посчитать вот такое выражение:

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

Дело в том, что наша система обнаружила ошибку, но продолжила считать, выдавая на выход блока последнее значение, которое было верным, а это 0. В принципе, предупреждение об ошибке система выводит, что хорошо. А еще лучше, что на это предупреждение можно не обращать внимания в данной конкретной задаче.

Первый вариант решения готов! Задачу мы решили, но не совсем честно, дело в том, что мы воспользовались тем, что в условиях задачи угол был 180 и она была симметрична, если же смотреть на условие задачи в общем, то кроме двух окружностей нам нужно было бы строить еще линии сечения и интегрировать уже не две простых окружности, а рассчитывать точки пересечения и подключать интеграторы. Как, например, на рисунке 13.

В этом случае сначала нужно интегрировать от a до b по прямой линии, а от b до с по линии окружности. Это можно, но схема потеряла бы свою простоту наглядность и понятность. Поэтому мы пойдем другим путем, начнем стрелять.

Ваше слово, товарищ Маузер

Прежде чем начать стрелять, давайте нарисуем мишень. У нас в задаче радиус внешней окружности 1 метр. Логично в качестве мишени использовать квадрат со стороной два метра, в которую впишется круг радиусом 1 метр.

Стрелять мы будем с помощью двух блоков «Равномерный шум». Один из них будет на каждом шаге выдавать новую случайную координату Х, второй — случайную координату Y в диапазоне от -1 до 1. Результаты работы этих блоков мы будем записывать в сигналы проекта с именем Хp и Yp.

Для того чтобы наглядно видеть, куда мы попали, сделаем прямо на схеме с использованием примитивов квадрат, который и будет нашей мишенью. Для этого я произвольным образом, используя четыре линии из панели блоков «Примитивы», нарисовал квадрат прямо на рабочем поле.

Для того чтобы убедиться, что фигура — это ровно квадрат, у каждой линии я настроил в свойствах координаты точек, у меня получился квадрат со стороной 250 условных единиц размеров рабочего поля.

Эти координаты позволят нам пересчитать Xp, Yp выстрелов таким образом, чтобы они все отображались в нарисованной мишени в правильных позициях.

Визуальное программирование — это хорошо, но скрипт лучше!

Поэтому перейдем к программированию в скрипте SimInTech, будем делать то, что ChatGPT еще не умеет делать, не знает она такого языка пока. Поэтому у нас пока есть счастливая возможность программировать сами, пользуемся, пока искусственный интеллект не отобрал эту возможность. Сам скрипт приведен на рисунке 16.

В окне скрипта заданы константами значения координат центра квадрата, нарисованного на схеме, и радиус вписанной окружности.

Эти значения позволяют нам обеспечить масштабирование координат, которые мы получаем на схеме в диапазоне от -1 до 1. Мы их просчитываем таким образом, чтобы они были отображены в мишени. Сама мишень — это квадрат 250 на 250 условных экранных единиц, а значит радиус круга, который можно вписать, будет 125 условных экранных единиц. Данный расчет выполняет функция ScalePoint (Relative: point).

Возможно, посмотрев на эту функцию, учитель математики опять будет в шоке, откуда при пересчете координат на плоскости появились функции комплексного числа real и imag?

На самом деле все просто, для координат точки мы используем тип переменных point, который содержит два числа. Это тип ровно такой же, какой мы используем и для комплексного числа, где координата Х точки соответствует реальной части комплексного числа, а Y — мнимой.

Также надо заметить, что направления осей на экране у нас соответствуют направлениям осей, принятым в Windows. А не так, как привыкли православные. Ось Y в богомерзкой Windows направлена вниз, поэтому при пересчете координат ее появляется знак минус.

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

За это отвечает функция function CreateCirclePoint(CCenter:point, R:integer), она создает кружок с заданным радиусом зеленого цвета.

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

Для того чтобы скрипт брал данные со схемы, нам необходимо добавить в список сигналов новые переменные: координаты точек Xp, Yp и количество точек PointCount. В меню Сервис главного окна есть пункт Сигналы. Войдите в диалоговое окно и добавьте эти сигналы, как показано на рисунке 18.

В настройках проекта установите конечное время расчета 5 и минимальный и максимальный шаг равными 1 секунде (см. рисунок 19).

Запускаем на расчет, и при каждом запуске программы в мишени появляются отметки от выстрелов. Координаты, выдаваемые блоками белого шума, записываются в сигналы проекта, на основании этих сигналов на каждом шаге проект скрип осуществляет отрисовку попаданий в мишени.

Отстрелялся - убери за собой!

Сейчас каждый новый расчет добавляет кружки к существующим, поле мишени зарастает зелеными кружочками, как пруд ряской. И каждый раз после расчета нужно удалять эти кружки вручную. А это не наш метод, пусть SimInTech чистит за нами поле.

Добавим в скрипт секцию инициализации. Поскольку при добавлении на схему примитивов их имена формируются путем добавления единицы к базовому имени, мы можем легко по имени их найти и удалить. Также обнуляем количество точек для нового расчета (см. рис. 21).

Если несколько раз запустить на расчет модель, видно, что один кружочек не меняет своего местоположения, это связано с тем, что его имя FillCircle и он не находится в цикле поиска, там нет единицы. Воспользуемся этим, вынесем его из мишени и сменим ему цвет на синий, этим кружочком мы будем отмечать наш центр масс, он будет у нас постоянно на схеме. После этого можно поставить время расчета 10 и шаг 0.1 и получить сотню отметок на мишени. При каждом запуске они будут менять свое положение.

Для любого действия нужна цель

Сейчас все выстрелы попадают в поле бесцельно. Давайте нарисуем в поле цель. Например, пусть считается попаданием те выстрелы, которые попадают в круг радиусом 1 м. Для этого мы, используя координаты точки, выполним расчет радиуса, на котором она находится от центра, и если радиус меньше заданного, то установим флаг попадания в цель. Схема расчета будет примерно, как на рисунке 22.

Добавим в список сигналов флаг HitTarget и будем его использовать в скрипте для того, чтобы поменять цвет у отметки, если она попала в круг. Добавим также счетчик попаданий в сигналы HitCount и будем его увеличивать при попадании. Не забудем также обнулить его при инициализации.

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

Если теперь увеличить время обстрела до 100 секунд, то мы получим 1000 выстрелов, а если до 1000 секунд, то будет 10000 выстрелов, мишень будет выглядеть как показано на рисунке 25.

После того как фигура цели наглядно поражена, можно начать считать площадь поражения. Вы уже догадались, как мы посчитаем площадь фигуры? Ведь у нас уже всё для этого есть!

Логика расчета простая: мы знаем, что обстрел ведется квадрата со стороной 2 м, тогда площадь обстрела составляет 4 м2. Поскольку обстрел ведется случайно и равномерно, количество попаданий в цель будет пропорционально площади цели. Посчитав соотношение попаданий к общим выстрелам, мы получим площадь цели.

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

И, как мы видим, всё прекрасно работает, площадь круга при R = 1 равна числу, и именно это число мы видим в результате расчета, с заданной точностью, конечно (см. рис. 26).

Цель поражена, пора делить шкуру

Теперь, формируя условия поражения цели по X и Y, мы можем легко формировать необходимую нам форму. Например, добавив в условия поражения радиус внутреннего круга, мы получим кольцо вместо круга.

Если для внутреннего круга r = 0.9, то схема и результат расчета будут как на рисунке 27.

Чтобы отрезать нижнюю половину и получить фигуру, как в исходной задаче, нужно добавить условие Y>0 (см. рис. 28).

Кстати, поскольку у нас тут дуэль структурных моделей, то, глядя на схему, я вижу, что мы для X и для Y делаем одинаковые математические действия, а это значит, что мы можем использовать векторность наших блоков. И вместо отдельных линий расчета х² и х² сделать всё в одной, но векторной. И тогда наша схема будет выглядеть вот примерно, как на рисунке 29. Источник белого шума выдает два числа Х и Y.

Перейдем теперь к расчету центра масс полученной фигуры. Используем тот же самый принцип подсчета точек, но теперь будем учитывать не только их количество, но и расстояние до центра. В самом деле, если считать, что все точки одинаково весят, то их центр масс будет где-то между ними, а это «между ними» легко посчитать как среднее арифметическое их координат.

Если считать в скрипте, то достаточно добавить два сигнала в список сигналов и две строчки в скрипт расчета.

Запускаем на расчет и видим, что после расчета у нас в списке сигналов проекта находятся координаты центра масс нашей фигуры. И, как мы видим, результат достаточно близкий к ответу, полученному по формуле. Х_с на миллиметр не дотягивает до 0, а на 3 мм меньше заданного (см. рис. 31).

Чтобы украсить задачу, вспомним про кружок центра масс, который мы сохранили в начале, и поместим его в рассчитанный центр масс на мишени. Для этого добавим скрипт перемещения фигуры, если мы попали в цель (см. рисунок 32).

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

Точность - вежливость королей и снайперов

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

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

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

Вывод простой, надо больше стрелять. Провел пять расчетов с 10000 выстрелов, получил такие вот пары количества попаданий и центра положение центра масс:721 - 0.612; 724 - 0.614; 734 - 0.586; 763 - 0.597; 810 - 0.612; 748 - 0.593.

Если стрелять 1 000 000 раз то у меня получилось 0.0606

Из всех искусств для нас важнейшим является рисование

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

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

Разберемся, как она работает. Для расчета нам нужна координата текущего выстрела X, Y и результат проверки на попадание в цель HitTarget. Можно было бы провести линию связи вокруг, но я использовал специальный блок переноса линий через переменную Hit. Мне так кажется аккуратнее. Пойдем последовательно.

1)   После определения факта попадания (0 — промазали, 1 — попали в цель), мы его перемножаем с координатой выстрела. Перемножая координату с результатом, мы либо обнуляем ее, выстрел мимо, и учитывать эту координату в среднем не нужно, либо не изменяем, поскольку умножаем на 1.

2) После того как мы получили либо два 0, либо 1 и координату попадания, мы собираем из них вектор, чтобы они ходили вместе.

3) Далее идет сложение с новой пары с суммой, которую мы запомнили. При этом у нас отдельно складываются координаты, полученные из блока выстрела, и количество попаданий.

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

5) Далее мы разделим вектор на две части: сумму координат попаданий и их количество.

6) Делим сумму на количество и получаем среднее арифметическое.

И делаем мы это сразу по двум координатам Х и Y, в итоге получается расчет центра масс на плоскости. Теперь вся задача решена честно в графическом виде. 

Нет предела совершенству

Поскольку у нас дуэль на лучшую схему, я решил еще упростить модель для решения данной задачи. Как было показано выше, при нахождении центра масс выстрелы мимо цели вообще не влияют, так зачем же мы пуляем в нижнюю часть площади мишени? Если можно прямо в настройках блока белого шума указать диапазон, и у нас исчезает один блок проверки для детекции попадания в цель. И вспомним, что у нас задача симметричная и нам не надо находить центр массы по Х, а значит, и здесь можно сократить диапазон обстрела.

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

Итоговая модель поиска центра масс методом массового обстрела по площади выглядит, как показано на рисунке 36. Я уменьшил радиус точек и не отображаю точки, которые не попадают в цель.

И результат получается значительно быстрее, процесс поиска отображен на графике, видно, что к приемлемой точности мы пришли где-то уже после 200 секунд расчета, это 2000 выстрелов (см. рис. 37).

Если закоментировать отображение в поражений в скрипте, то можно зарядить и миллион выстрелов, у меня при этом получилось - 0.60532902

Модель расчета можно взять здесь ...

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

Как я говорил в начале, у нас здесь идет дуэль, поэтому выкладываю схемы расчета центра масс в виде структурной схемы от Валерия Федоровича и выношу на голосование: какая схема красивее, четыре варианта решения задачи в виде струкутрной схемы: