Что можно улучшить в устройстве из 60-х?

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

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

Источник: http://www.prom59.ru/sushilka.htm

Если сделать срез, то перед нами предстанет такая картина

Если описывать простыми словами, то в барабане мокрый песочек и горячий воздух, барабан вращается, лопатки внутри загребают песочек и несут его вверх, откуда песочек ссыпается вниз. Результатом этого процесса является сухой песочек.
Устройства эти применяются уже давно и чтобы чем-то удивить производителей, нужно сильно постараться. Но мы попробуем.
Основной задачей таких машин является сушка и чем быстрее она будет проходить, тем лучше. Чтобы улучшить производительность нужно что-то поменять в этой машине. Но барабан - он и в Африке барабан, менять можно только лопатки.
Решено было сделать простую симуляцию системы частиц во вращающемся барабане. Частиц очень много, поэтому симуляцию хорошо было бы производить на CUDA, тем более ноут с нужной видеокарточкой давно ждал своего часа.
Качаем драйвера и SDK, ставим. В примерах видим подходящий шаблон, который называется particles. Вот на youtube нашел видео с записью процесса симуляции

Здесь, как мы видим, все происходит в неподвижном кубе. Нам нужно переделать куб на цилиндр с конфигурируемыми лопатками, а потом заставить это вращаться.
В проекте открываем файл particles_kernel.cu и ищем там описание структуры integrate_functor. Здесь содержится описание процесса интегрирования. А процесс интегрирования включает в себя не только определение новой координаты через значение силы, но и проверку граничных условий.

if (pos.z > 0.25f - params.particleRadius) { pos.z = 0.25f - params.particleRadius; vel.z *= -0.6f; }

Здесь мы видим обработку события столкновения частицы с одной из стенок куба. Это условие очевидно. Если частица в результате изменения координаты оказалась за границей стенки, то мы ее перемещаем так, чтобы она находилась ровно на стенке, а одну из составляющих ее скорости меняем на противоположную и учитываем потерю энергии при соударении.
При замене куба на цилиндр эта стенка как раз останется, поэтому здесь мы менять ничего не будем, не будем менять и обработку соударения с противоположной стенкой. А вот другие обработки закомментируем.
Столкновение с цилиндром будет происходить, когда проекция радиус-вектора частицы на плоскость XY будет по модулю больше единицы (цилиндр единичного радиуса, и его ось лежит на оси Z).
Для такой проекции мы в начале функции объявляем переменную post

float3 post = make_float3(posData.x, posData.y, 0.0f);

А потом проверяем

if (length(post) > 1.0f) { 

Что же нужно делать в случае соударения частицы с внутренней поверхностью цилиндра? Будем раскладывать скорость по двум осям: одна ось образована радиусом цилиндра, а другая перпендикулярна ей и сонаправлена с вектором тангенциальной скорости. По нормальной оси будет происходить то же, что и при соударении со стенкой (составляющая скорости вдоль радиус-вектора меняет напрвление, а модуль уменьшается из-за частичного перехода кинетической энергии в тепловую), а по тангенциальной будет действовать сила трения, зависящая от силы реакции опоры.

if (length(post) > 1.0f) {
    post = post / length(post); //возвращаем частицу на стенку цилиндра
    pos.x = post.x;
    pos.y = post.y;
    norm_axis = post; //нормальная ось есть радиус-вектор
//тангенциальная ось вычисляется исходя из предположения, что
//срость не будет направлена только вдоль радиус-вектора
//всегда должна быть хоть какая-то составляющая,
//которая и равна vel - dot(vel, norm_axis)*norm_axis
    tan_axis = (vel - dot(vel, norm_axis)*norm_axis) /
               length(vel - dot(vel, norm_axis)*norm_axis);
    float tan_vel = dot(vel, tan_axis);
    force = -1.5f*dot(vel, norm_axis)*norm_axis/deltaTime;
    force -= friction_coef*length(force)*tan_axis;
    tan_axis.x = norm_axis.y;
    tan_axis.y = -norm_axis.x;
    tan_axis.z = 0.0f;
}

К этому моменту уже нужно провести эксперимент: проверить, "работает" ли сила трения. Для этого меняем коэффициент трения и смотрим, как "устаканятся" частицы.

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

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

Если вы еще не поняли, то лопатка в нашей модели представляется как набор четырехугольников (если на них смотреть сбоку). На картинке изображен сегмент, прилегающий к стенке барабана, и часть следующего сегмента. Толщина лопатки такова, что за один шаг симуляции никакая из частиц не пролетит и половины толщины лопатки (ну это если ко всему подходить систематично и серьезно, на самом деле выбираем так, чтобы выглядело красиво и этого окажется более чем достаточно).
Узнать, столкнулась ли частица с лопаткой, можно, проверив, находится ли она внутри прямоугольника сегмента. Если да, то обрабатываем столкновение.
В результате столкновения частицу нужно будет перенести к одной из стенок. Чтобы узнать к какой стенке переносить частицу, смотрится ее предыдущее положение. Для этого из радиус-вектора текущего положения вычитается вектор скорости, умноженный на шаг времени. А потом проверяется, по разные стороны от какой плоскости лопатки находятся эти две точки. И к той плоскости и переносим частицу. Частицы переносятся по перпендикуляру к плоскости. Есть одна тонкость для сегмента, изображенного на рисунке выше. В некоторых случаях частицы начнут переноситься за границы цилиндра. И мы будем наблюдать "залипания" частиц. Поэтому здесь необходим дополнительный перенос вдоль направляющей сегмента до тех пор, пока частица не окажется на поверхности цилиндра.
Вот код, выполняемый при столкновении с первым сегментом.

if ( InsideTriangle(post, params.p1[i], params.p2[i], params.p3[i]) ||
     InsideTriangle(post, params.p2[i], params.p3[i], params.p4[i]) ) {
    prev_pos = pos - vel * deltaTime*4;
    if (DiffSides(pos, prev_pos, params.p1[i], params.p3[i])) {
        D = -(params.A1[i]*pos.x+params.B1[i]*pos.y+params.C11[i])/
             (params.A1[i]*params.A1[i]+params.B1[i]*params.B1[i]);
        pos.x += D*params.A1[i];
        pos.y += D*params.B1[i];
        norm_axis = make_float3(params.A1[i], params.B1[i], 0.0f);
        norm_axis = norm_axis/length(norm_axis);
        vel -= 1.5f * dot(norm_axis, vel)*norm_axis;
        post = pos;
        post.z = 0.0f;
        if (length(post) > 1.0f) {
            tan_axis = make_float3(-params.B1[i], params.A1[i], 0.0f);
            D = dot(tan_axis, post)*dot(tan_axis, post) -
                length(post)*length(post) + 1.0f;
            float ksi1 = (-dot(post, tan_axis) + sqrt(D));
            float ksi2 = (-dot(post, tan_axis) - sqrt(D));
            if(abs(ksi1) < abs(ksi2)) {
                pos += ksi1 * tan_axis;
            }
            else {
                pos += ksi2 * tan_axis;
            }
        }
}

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

Вот много односегментных лопаток.

Вот четырехсегментные.

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

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

Ну что сказать... Похоже.
Ну и видео еще вставим.

Код проекта можно будет найти в конце, а пока перейдем к оптимизации процесса.
Оптимизация - это довольно сложная область математики, имеющая в своем арсенале кучу методов. Критерии выбора этих методов известны профессионалам. Поэтому сразу стоит сказать, что наш метод (эволюционная оптимизация) может и не являться оптимальным. Выбран он был, потому что было просто интересно посмотреть, что из этого получится.
В Интернете много статей по применению метода.

Вкратце:

1. Есть популяция лопаток (пусть будет их 20 штук) - левый верхний угол
2. Из этой популяции выбираются пары для скрещивания - получается новая популяция скрещенных лопаток.
3. На некоторые параметры получившихся лопаток может действовать мутация (с низкой вероятностью).
4. Из множества имеющихся лопаток выбираются лучшие - те, которые дают наивысшее значение целевой функции.
5. Процесс повторяется.

С точки зрения математики происходит поиск экстремумов целевой функции. В нашем случае целевая функция является функцией 10 параметров (5 длин сегментов и 5 углов между ними). Так что поиск происходит в 10-мерном пространстве. Часто бывает, что все точки стремятся к одному экстремуму. Если нужно найти хотя бы один, то можно дальше не выпендриваться. А если нужно искать несколько, то применяют мутации, которые "выбивают" точку из окрестности экстремума и в результате чего поиск может перекинуться на другой экстремум.

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

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

struct is_inside
{
    __host__ __device__
    bool operator() (const float4& x)
    {
        float2 d_pos3 = make_float2(x.x, x.y);
        return (length(d_pos3) < params.min_rad)&&
               (d_pos3.y>0.15)&&(d_pos3.y>-d_pos3.x/2+0.15);
    }
};

int count_particles(float* pos, uint num_part, float rad)
{
    thrust::device_ptr<float4> d_pos4((float4 *)pos);
    return thrust::count_if(d_pos4, d_pos4+num_part, is_inside());
}

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

Как видно, количество частиц в конце концов выравнивается. Вот поподробнее в установившемся состоянии.

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

def CreateFirstGeneration(generation, results, generation_size) :
    while len(results) < generation_size :
        pers = GenerateRandomPerson()
        if TestAppend(pers, generation, results) == 0 :
            print(str(len(results)))
        print("Finally first gen:\n" + str(results))

def SaveGeneration(generation, results, generation_count) :
        genf = open("history/gener" + str(generation_count) + ".txt", "w+")
        for p in generation :
                for g in p[0] :
                        genf.write(str(g))
                        genf.write("\n")
                for g in p[1] :
                        genf.write(str(g))
                        if p != generation[-1]:
                                genf.write("\n")
                        elif g != p[1][-1] :
                                genf.write("\n")
                if p != generation[-1] :
                        genf.write("\n")
        genf.close()
        
        resf = open("history/results" + str(generation_count) + ".txt", "w+")
        for r in results :
                resf.write(str(r))
                if r != results[-1] :
                        resf.write("\n")
        resf.close()

def CreateNewGeneration(_prev_gen, _prev_res, _new_gen, _new_res) :
        print("Creating new generation")
        SortPersons(_prev_gen, _prev_res)
        while len(_new_gen) < len(_prev_gen) :
                father = randint(0, len(_prev_gen)-1)
                mother = randint(0, len(_prev_gen)-1)
                print("Creating child" + str(father) + "+" + str(mother))
                child = GenerateChild(_prev_gen[mother], _prev_gen[father])
                TestAppend(child, _new_gen, _new_res)
        SortPersons(_new_gen, _new_res)
        _new_gen[:] = _prev_gen[0:int(len(_prev_gen)/2)]+_new_gen[0:int(len(_new_gen)/2)]
        _new_res[:] = _prev_res[0:int(len(_prev_res)/2)]+_new_res[0:int(len(_new_res)/2)]
        print("Finaly:\n" + str(new_res))

def GenerateChild(pers1, pers2) :
        child = []
        child.append([])
        child.append([])
        for i, gene in enumerate(pers1[0]) :
                alpha = uniform(-0.25, 1.25)
                child[0].append(gene + alpha*(pers2[0][i]-gene))
        for i, gene in enumerate(pers1[1]) :
                alpha = uniform(-0.25, 1.25)
                child[1].append(gene + alpha*(pers2[1][i]-gene))
        return child

def TestAppend(person, generation, results) :
        if CanItLive(person) :
                res = RunSimulation(person)
                if res > 20.0 :
                        generation.append(person)
                        results.append(res)
                        print(str(results))
                        return 0
                else :
                        print("Person can't live : it works bad")
        else :
                print("Person can't live : it's cripple")
        return 1

def RunSimulation(person) :
        WritePersonToFile(person)
        LaunchApp()
        return GetParticlesNumber()

def GetParticlesNumber() :
        num = [int(line.strip()) for line in open("result.txt")]
        result = sum(num[int(-len(num)/5*2):-1])/(int(len(num)/5*2))
        return result

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

cur_generation = []
prev_generation = []
cur_results = []
prev_results = []
generation_counter = 5
#print("Creating first generation")
#CreateFirstGeneration(prev_generation, prev_results, 40)
ResumeProcessing(prev_generation, prev_results, generation_counter)
generation_counter += 1
while generation_counter < 10 :
        SaveGeneration(prev_generation, prev_results, generation_counter)
        CreateNewGeneration(prev_generation, prev_results, cur_generation, cur_results)
        prev_generation = cur_generation
        cur_generation = []
        prev_results = cur_results
        cur_results = []
        generation_counter += 1
SaveGeneration(prev_generation, prev_results, generation_counter)

Так как все для домашнего использования, то я не стесняюсь комментировать/раскомментировать некоторые строки в зависимости от ситуации. Так можно возобновить поиск с предыдущего поколения, если машина перегрелась или мамка выгнала спать надоело наблюдать за процессом.
Качайте вот это Polygon.zip (1022 KB), ставьте интерпретатор Python, запускайте, меняйте количество особей в поколении и наблюдайте за чудесами. В архиве прилагается исполняемый файл, так что с первой частью можно не разбираться, если нет желания. И да, прогал под Windows, если хотите под Linux, то нужно поменять директории (я где-то использовал абсолютные адреса).
Какой же результат? Оказался вполне ожидаемым. Уже лет 30 назад наинженерили без всяких моделирований оптимальные формы лопаток.

В нашем случае получилось то же. Только вот первый вариант вышел даже слегка лучше. Второй сегмент скруглился, тем самым увеличив количество единовременно поднимаемых частиц. Бессознательный поиск дал результат, который несет в себе довольно "глубокую" мысль!
Еще одна интересная деталь бессознательного поиска. Длины сегментов и кривизна лопаток оказываются в процессе эволюции такими, что "кучки" частиц, которые собираются на лопатке оказываются по размеру ровно такими, что ссыпаться они начинают едва касаясь следующей лопатки. Замечательно, я считаю!
Ну напоследок еще картинок.

Обещанные исходники симулятора. particles.zip (13.4 MB) Чтобы пробовать разные лопатки можете подредактировать, как вам хочется, файл params.txt.
В файле конфигурации находятся сначала углы между продолжением сегмента и следующим сегментом, а потом длины сегментов, каждый на новой строке.
И простите за не очень удачный код, который я поленился в итоге переделывать на switch, который выглядел в итоге так

Не делайте так.