Autodif types in PyGMO optimisation problems¶

Introduction¶

This example shows a basic usage of autodif types in optimisation problems with PyGMO and a comparison with numerical derivatives.

In [1]:

import pygmo as pg

Import GODOT automatic differentiation module.

In [2]:

from godot.core import autodif
from godot.core.autodif import bridge as br

Import other modules.

In [3]:

import math
import numpy as np
import time

PyGMO problem class using gradient from numerical differences¶

This example problem is taken from the PyGMO online documentation and is detailed here: https://esa.github.io/pygmo2/tutorials/coding_udp_constrained.html

In [4]:

class my_udp:
    def fitness(self, x):
        obj = 0
        for i in range(3):
            obj += (x[2*i-2]-3)**2 / 1000. - (x[2*i-2]-x[2*i-1]) + math.exp(20.*(x[2*i - 2]-x[2*i-1]))
        ce1 = 4*(x[0]-x[1])**2+x[1]-x[2]**2+x[2]-x[3]**2
        ce2 = 8*x[1]*(x[1]**2-x[0])-2*(1-x[1])+4*(x[1]-x[2])**2+x[0]**2+x[2]-x[3]**2+x[3]-x[4]**2
        ce3 = 8*x[2]*(x[2]**2-x[1])-2*(1-x[2])+4*(x[2]-x[3])**2+x[1]**2-x[0]+x[3]-x[4]**2+x[0]**2+x[4]-x[5]**2
        ce4 = 8*x[3]*(x[3]**2-x[2])-2*(1-x[3])+4*(x[3]-x[4])**2+x[2]**2-x[1]+x[4]-x[5]**2+x[1]**2+x[5]-x[0]
        ci1 = 8*x[4]*(x[4]**2-x[3])-2*(1-x[4])+4*(x[4]-x[5])**2+x[3]**2-x[2]+x[5]+x[2]**2-x[1]
        ci2 = -(8*x[5] * (x[5]**2-x[4])-2*(1-x[5]) +x[4]**2-x[3]+x[3]**2 - x[4])
        return [obj, ce1,ce2,ce3,ce4,ci1,ci2]

        return np.asarray(f)

    def get_bounds(self):
        return ([-5]*6,[5]*6)
    def get_nic(self):
        return 2
    def get_nec(self):
        return 4
    def gradient(self, x):
        return pg.estimate_gradient_h(lambda x: self.fitness(x), x)

Create user-defined problem instance for PyGMO.

In [5]:

prob = pg.problem(my_udp())

Setup NLOPT optimisation algorithm for PyGMO.

In [6]:

algo = pg.algorithm(uda = pg.nlopt('auglag'))
algo.extract(pg.nlopt).local_optimizer = pg.nlopt('var2')
algo.set_verbosity(200)                    # in this case this correspond to logs each 200 objevals

Setup population with one individual and initial guess.

In [7]:

pop = pg.population(prob, size=0)  # empty population
pop.push_back( [0,0,0,0,0,0] )     # add initial guess
pop.problem.c_tol = [1E-6] * 6     # set the constraints tolerance

Run the optimisation algorithm (i.e. "evolve the population" in PyGMO jargon) and print the results.

In [8]:

start = time.time()
pop = algo.evolve(pop)    # run the optimisation
print("Execution time [sec] = ", time.time() - start)
print("Solution decision vector = ", pop.champion_x)
print("Solution cost function = ", pop.champion_f[0])
print("Solution constraints = ", pop.champion_f[1:])

Execution time [sec] =  4.616239786148071
Solution decision vector =  [ 0.46422913  0.45459031 -0.23568187  0.40464072  0.31477216  0.70663412]
Solution cost function =  2.2596649939249227
Solution constraints =  [ 1.45909206e-08 -1.03483065e-09  6.16791684e-09 -2.64623740e-09
 -8.18679014e-01 -9.11250519e-11]

 objevals:        objval:      violated:    viol. norm:

         1          3.027              4         5.4641 i
       201        1.14765              5       0.355196 i
       401        1.47742              5      0.0964976 i
       601        2.02658              4      0.0190178 i
       801        2.12823              5       0.011737 i
      1001        2.23834              5     0.00161745 i
      1201        2.25766              5    0.000162072 i
      1401        2.38666              5       0.031823 i

Optimisation return status: NLOPT_XTOL_REACHED (value = 4, Optimization stopped because xtol_rel or xtol_abs was reached)

PyGMO Problem class using autodif types¶

The transcribed problem class replaces the arithmetic operations with functions from the autodif.bridge module for type-agnostic operations. The evaluate() function can then return both numpy arrays or autodif vector containing derivatives w.r.t. the decision vector, x. The gradient() method extracts those derivatives and returns them in a format as expected by PyGMO.

In [9]:

class my_autodif_udp:
    def fitness(self, x):
        out = self.evaluate( x, False )
        return out

    def gradient( self, x ) :
        out = self.evaluate( x, True ) 
        grad = out.gradient()[ out.parameters()[0] ] # extract gradient w.r.t. parameters from autodif.Vector
        return grad.flatten( order = 'C' )           # reshape as expected by PyGMO

    def evaluate( self, x_float, withPartials ) :

        if withPartials:
            x = autodif.Vector(x_float, "x")  # autodif vector type with x as "leaf"
        else:
            x = x_float

        obj = 0
        for i in range(3):
            ind = np.mod(2*i-2, 6)  # needed since currently autodif.Vectors can only be addressed with positive indices
            obj +=  br.pow( x[ind]-3, 2) / 1000. - (x[ind]-x[ind+1]) + br.exp(20.*(x[ind]-x[ind + 1]))
        ce1 = 4*br.pow(x[0]-x[1],2)+x[1]-br.pow(x[2],2)+x[2]-br.pow(x[3],2)
        ce2 = 8*x[1]*(br.pow(x[1],2)-x[0])-2*(1-x[1])+4*br.pow(x[1]-x[2],2)+br.pow(x[0],2)+x[2]-br.pow(x[3],2)+x[3]-br.pow(x[4],2)
        ce3 = 8*x[2]*(br.pow(x[2],2)-x[1])-2*(1-x[2])+4*br.pow(x[2]-x[3],2)+br.pow(x[1],2)-x[0]+x[3]-br.pow(x[4],2)+br.pow(x[0],2)+x[4]-br.pow(x[5],2)
        ce4 = 8*x[3]*(br.pow(x[3],2)-x[2])-2*(1-x[3])+4*br.pow(x[3]-x[4],2)+br.pow(x[2],2)-x[1]+x[4]-br.pow(x[5],2)+br.pow(x[1],2)+x[5]-x[0]
        ci1 = 8*x[4]*(br.pow(x[4],2)-x[3])-2*(1-x[4])+4*br.pow(x[4]-x[5],2)+br.pow(x[3],2)-x[2]+x[5]+br.pow(x[2],2)-x[1]
        ci2 = -(8*x[5] * (br.pow(x[5],2)-x[4])-2*(1-x[5]) +br.pow(x[4],2)-x[3]+br.pow(x[3],2) - x[4])
        return br.concatenate([obj, ce1,ce2,ce3,ce4,ci1,ci2])

    def get_bounds(self):
        return ([-5]*6,[5]*6)
    def get_nic(self):
        return 2
    def get_nec(self):
        return 4

Create user-defined problem instance with autodif types.

In [10]:

prob = pg.problem(my_autodif_udp())

Re-set initial decision vector and evolve the population.

In [11]:

pop = pg.population(prob, size=0)  # empty population
pop.push_back( [0,0,0,0,0,0] )     # add initial guess
pop.problem.c_tol = [1E-6] * 6     # set the constraints tolerance

start = time.time()
pop = algo.evolve(pop)    # run the optimisation
print("Execution time [sec] = ", time.time() - start)
print("Solution decision vector = ", pop.champion_x)
print("Solution cost function = ", pop.champion_f[0])
print("Solution constraints = ", pop.champion_f[1:])

Execution time [sec] =  1.1535863876342773
Solution decision vector =  [ 0.46422915  0.45459032 -0.23568187  0.40464072  0.31477224  0.70663416]
Solution cost function =  2.2596653630417762
Solution constraints =  [ 1.62131411e-08 -7.00207578e-09  2.76815351e-08 -1.71066832e-08
 -8.18678996e-01 -6.62564559e-10]

 objevals:        objval:      violated:    viol. norm:

         1          3.027              4         5.4641 i
       201        1.14768              5       0.355326 i
       401        1.44876              5       0.118643 i
       601        2.02301              5      0.0198603 i
       801        2.18888              5      0.0054834 i
      1001        2.23798              5     0.00163888 i
      1201        2.25318              5    0.000486845 i
      1401        2.25907              5    4.20613e-05 i
      1601        2.25966              0              0

Optimisation return status: NLOPT_XTOL_REACHED (value = 4, Optimization stopped because xtol_rel or xtol_abs was reached)

The same local minimum has been reached. The execution time has gone down from 7.64 sec to 2.51 sec. But more importantly, the autodif derivatives are more reliable than the numerical ones where the user has to choose the differencing step parameter. Especially, for real-world problems where optimisation parameters of different scales are involved, numerical differences become more tricky to deal with.