ТАУ-Дарвинизм: реализация на Ruby

Предисловие

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

Мультфильм «Пластилиновая ворона».

В этой статье представляю Вашему вниманию реализацию эволюционного подхода к идентификации динамической системы. Программа написана на языке Ruby версии 1.9.2 (gems: NArray, GNUPlot). Заглянув под кат найдете пример вещественного кодирования генной информации и подходящего для него алгоритма скрещивания («flat crossover»).

Введение

Начну с описания задачи «черного ящика». Предположим у нас есть некий объект, передаточную функцию которого мы бы хотели узнать (упрощенную, близкую или эквивалентную к реальной). Все, что мы можем сделать, это подать на вход объекту некоторое воздействие и измерить некоторый выходной параметр (например, выходное напряжение). Как нам подглядеть, что у ящика внутри? Иными словами: как определить примерную математическую модель исследуемого объекта?
В теории автоматического управления (ТАУ) данная задача называется идентификацией. Существуют разные способы идентификации, основанные на серьезной математике и сложные в реализации. Однако есть возможность получить приближенную математическую модель сравнительно простым способом, а именно используя стохастический подход с применение эволюционных алгоритмов.
Подробности Вы можете прочитать в первой статье, посвященной данному вопросу: ТАУ-Дарвинизм, а здесь я остановлюсь на программной реализации.

Реализация

При написании алгоритмов использована библиотека «gga4r» (rubygems.org). Данная библиотека была переписана для приведение в работоспособное состояние для совместимости с ruby 1.9. Автор этой библиотеки (спасибо ему большое!) реализовал подход эволюционной обертки над сущностями, реализующими представления различных по содержанию, но одинаковых по форме эволюционирующих особей. Предлагается реализовывать свой собственный класс, имплементирующий генетические операторы. У меня получилось это так:

class TFIndivid
    attr_reader :filter, :fitness, :dna, :num_len, :denom_len
    def initialize(num,denom,k=1.0,ts=@@ts)
      @num_len,@denom_len,@filter = num.length,denom.length,SSF::tf2ssd(num, denom, k, ts)
      @fitness = 1/(step_response_deviation(@filter)+@@eps)
      @fitness = 0.0 if @fitness.infinite?
      @dna = num + denom + [k] + [ts]
      self
    end
    def integrity_ok?(dna)
      return dna.length == @num_len + @denom_len + 2
    end
    def TFIndivid.flat_crossover(par1, par2)
      r = Random.new
      child = []
      par1.size.times { |i|
        if par1[i]>par2[i]
          child << r.rand(par2[i]..par1[i])
        else
          child << r.rand(par1[i]..par2[i])
        end
      }
      child
    end
    def recombine(other)
      child = TFIndivid::flat_crossover(self.dna, other.dna)
      res = self.integrity_ok?(child) ? TFIndivid.new(child[0..@num_len-1], child[@num_len..-3], child[-2], child[-1]) : []
      [res]
    end
    def mutate
      rnd = Random.new
      mutate_point = rnd.rand(@dna.size-1)
      @dna[mutate_point] *= rnd.rand(1.0)
    end
end

В примере на сайте «gga4r» класс особи наследовался от Array. Я решил поменять немного логику и завел дополнительное свойство "@dna". Сделано это было для удобства и использования конструктора, в котором происходит создание цифрового фильтра для данных коэффициента усиления, числителя и знаменателя передаточной функции (в непрерывном времени), а также периода квантования по времени. Также в конструкторе переданные параметры сливаются в один одномерный массив для работы с ними посредством генетических алгоритмов. Фитнесс особи вычисляется сразу в конструкторе и хранится в соответствующем свойстве.
Функция «flat_crossover» реализует упомянутый в самом начале статьи вещественного скрещивания. Это, наверное, самая простая его реализация. Существует еще несколько, в которых детей реализуется несколько (здесь генерируется один потомок) и которые мне наверняка предстоит реализовать. В данном же случае порождаемое дитя формируется путем выбора для каждого гена случайного числа в интервале между значениями генов родителей.
Оператор мутации реализован в виде умножения значения гена на случайное число в диапазоне [0.0..1.0]. Таким образом некоторые из генов «притягиваются» к нулю. Выбранный оператор скрещивания формирует в основном особи со значениями генов в средней части допустимого диапазона значений. Выбранный способ мутирования играет роль противовеса.
Сама обертка эволюционного алгоритма работает как описано в предыдущей статье ТАУ-Дарвинизм: формируется популяция, проводится селекция, рекомбинация (скрещивание), мутация. Внутренности сможете посмотреть в исходниках программы (ссылка в конце статьи). Ее пришлось немного переписать, как написано выше, для возможности использования с ruby 1.9. Но на самом деле библиотека у меня отказалась работать и с версией 1.8 на kubuntu 10.10 (не выяснял почему). Как бы то ни было, можете сравнить с оригиналом «gga4r» — изменений немного. Кое что еще предстоит изменить.
Фитнесс-функция реализована как среднеквадратическое отклонение между выходным сигналом тестового фильтра и сгенерированного по генам особи:

def step_response_deviation(filter)
  x, test_y,y = step_response(filter)
  av_err = 0.0
  100.times {|i|
    err = test_y[i] - y[i]
    av_err += err*err
  }
  Math.sqrt(av_err / 99)
end
def step_response(filter)
  @@test_filter.reset
  filter.reset
  x = []; 100.times { |i| x << i*@@ts }
  test_y,y = [],[]
  100.times {
    test_y << @@test_filter.step(@@test_input)
    y << filter.step(@@test_input)
  }
  return x,test_y,y
end

Это немного отличается от постановки задачи («дана запись выхода и входа черного ящика; требуется...»), но так удобнее, а суть та же.

Результаты

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

• Имеется эталонная (тестовая) передаточная функция — числитель: [1.5, 1.0], знаменатель: [0.75, 0.35, 1.1]
• Коэффициент усиления тестового фильтра выбрал равным единице (коэффициент передаточной функции при построении дискретного фильтра к нему приводится в функции SSF::tf2ssd — модуль «ssf.rb» в исходниках)
• Период квантования — 0.1 сек., для выбранного тестового фильтра вполне достаточно малый период
• Размер популяции — 350 в начале и 500 максимально
• Число итераций — 50

Оказалось, что не смотря на заданный максимум размер популяции превышал заданное значение. После запуска получил следующий результат:
$ ruby1.9.1 main.rb
Iteration: 1
Population size: 387
best Individ: [ 1.552605923140199 1.8418549861379885 4.4154048910142265 5.556332886121151 5.320061971918478 1.6334223860651258 0.1 ]
best fitness: 2.6508867715774027
mean fitness: 1.1722382322077871
mean fitness recalc 1.1722382322077871
******************************
Iteration: 2
Population size: 421
best Individ: [ 0.37594352736713244 0.41988779230454 6.092873656120998 3.575618643442387 8.786369744495486 1.3186424599689655 0.1 ]
best fitness: 2.9780762386062767
mean fitness: 1.381001802845182
mean fitness recalc 1.381001802845182
******************************
Iteration: 3
Population size: 456
best Individ: [ 5.355702236617921 2.9021864520499134 3.979973441887163 1.1587534831116912 5.818934614234973 0.6934218845420625 0.1 ]
best fitness: 4.008150118088937
mean fitness: 1.6192001539010286
mean fitness recalc 1.6192001539010286
******************************
Iteration: 4
Population size: 501
best Individ: [ 6.981596537433821 7.0034200556235495 6.592795290155172 3.1179936580946577 7.236199584387684 1.4664747534746492 0.1 ]
best fitness: 4.629268673426835
mean fitness: 1.8244067682429481
mean fitness recalc 1.8244067682429481
******************************
Iteration: 5
Population size: 537
best Individ: [ 5.917012124258587 6.182284010197045 5.214809836688884 3.9333556599339055 8.349166334410114 1.544147649767568 0.1 ]
best fitness: 4.900537478511452
mean fitness: 2.024810983956715
mean fitness recalc 2.024810983956715
******************************
Iteration: 6
Population size: 551
best Individ: [ 5.917012124258587 6.182284010197045 5.214809836688884 3.9333556599339055 8.349166334410114 1.544147649767568 0.1 ]
best fitness: 4.900537478511452
mean fitness: 2.3489882909253152
mean fitness recalc 2.3489882909253152
******************************
Iteration: 7
Population size: 540
best Individ: [ 7.667944581517669 7.685145812466274 3.6039889820553768 2.4390348703480536 5.087422184655714 1.4689450977845646 0.1 ]
best fitness: 8.166839319133336
mean fitness: 2.652950357867813
mean fitness recalc 2.652950357867813
******************************
Iteration: 8
Population size: 552
best Individ: [ 7.667944581517669 7.685145812466274 3.6039889820553768 2.4390348703480536 5.087422184655714 1.4689450977845646 0.1 ]
best fitness: 8.166839319133336
mean fitness: 2.9467084852969725
mean fitness recalc 2.9467084852969725
******************************
Iteration: 9
Population size: 553
best Individ: [ 5.7806312313223485 4.681216450534634 5.384196754050039 3.0234885021495357 7.238140602489246 1.1815749271908358 0.1 ]
best fitness: 9.966075154360498
mean fitness: 3.335389533280396
mean fitness recalc 3.335389533280396
******************************
Iteration: 10
Population size: 543
best Individ: [ 5.7806312313223485 4.681216450534634 5.384196754050039 3.0234885021495357 7.238140602489246 1.1815749271908358 0.1 ]
best fitness: 9.966075154360498
mean fitness: 3.818330258524996
mean fitness recalc 3.818330258524996
******************************
Iteration: 11
Population size: 555
best Individ: [ 5.7806312313223485 4.681216450534634 5.384196754050039 3.0234885021495357 7.238140602489246 1.1815749271908358 0.1 ]
best fitness: 9.966075154360498
mean fitness: 4.241020711174743
mean fitness recalc 4.241020711174743
******************************
Iteration: 12
Population size: 563
best Individ: [ 7.114438733870453 5.7051128886012386 3.8495059720675098 1.8164317233081309 5.543227442940334 1.1352374472380862 0.1 ]
best fitness: 11.644557896655645
mean fitness: 4.980294253355436
mean fitness recalc 4.980294253355436
******************************
.......
******************************
Iteration: 23
Population size: 542
best Individ: [ 5.949029964618817 4.681516090696933 4.8007731550065325 2.562344672524173 6.969591955243219 1.1684490266387202 0.1 ]
best fitness: 20.468430876331446
mean fitness: 11.25191007878279
mean fitness recalc 11.25191007878279
******************************
......
******************************
Iteration: 36
Population size: 550
best Individ: [ 6.304563716617802 4.937028188768713 5.01636959411469 2.5019167234335344 7.143946485486306 1.1690443373254094 0.1 ]
best fitness: 30.361894930102864
mean fitness: 22.65217005076974
mean fitness recalc 22.65217005076974
******************************
........
******************************
Iteration: 46
Population size: 542
best Individ: [ 6.548169493382168 5.060607979367818 4.739775170940725 2.3444848651883947 6.741642980854007 1.1564890916088506 0.1 ]
best fitness: 32.826565607333514
mean fitness: 25.606090666108418
mean fitness recalc 25.606090666108418
******************************
Iteration: 47
Population size: 550
best Individ: [ 6.548169493382168 5.060607979367818 4.739775170940725 2.3444848651883947 6.741642980854007 1.1564890916088506 0.1 ]
best fitness: 32.826565607333514
mean fitness: 25.489911664932283
mean fitness recalc 25.489911664932283
******************************
Iteration: 48
Population size: 561
best Individ: [ 6.548169493382168 5.060607979367818 4.739775170940725 2.3444848651883947 6.741642980854007 1.1564890916088506 0.1 ]
best fitness: 32.826565607333514
mean fitness: 25.44926190395142
mean fitness recalc 25.44926190395142
******************************
Iteration: 49
Population size: 550
best Individ: [ 3.1439340250602816 5.060607979367818 4.739775170940725 2.3444848651883947 6.741642980854007 1.1564890916088506 0.1 ]
best fitness: 32.826565607333514
mean fitness: 26.108921277427932
mean fitness recalc 26.108921277427932
******************************
Iteration: 50
Population size: 559
best Individ: [ 3.1439340250602816 5.060607979367818 4.739775170940725 2.3444848651883947 6.741642980854007 1.1564890916088506 0.1 ]
best fitness: 32.826565607333514
mean fitness: 26.23635084771665
mean fitness recalc 26.23635084771665
******************************
Первые два числа в «best Individ» — числитель передаточной функции, следующие три — знаменатель, далее — масштабный коэффициент и период квантования. Наилучший фитнесс, равный 32 «попугаям» — достаточно хороший результат. В подтверждение следующий рисунок:

На нем изображены графики переходных процессов тестового и полученного в ходе эволюции фильтров при подаче им на вход единичного ступенчатого воздействия. «Кто» из графиков «Кто» пока не могу сказать — недопилил функциональность по работе с gnuplot, но это и не имеет значение в данном случае, т.к. невооруженным глазом видно, что результаты очень близки.

Заключение

Программа реализована в ходе исследовательских работ по программе «У.М.Н.И.К». Код сырой, но работающий (по крайней мере у меня). В планах оснастить программу неким простеньким GUI для удобства. Ниже ссылка на код.
UPD1: исходники размещены на GitHub.
Репозиторий на GitHub