Блог

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

Я планирую использовать ipopt для оптимизации позиции. В настоящее время я изучаю, как использовать ipopt, и именно здесь я сейчас сталкиваюсь с проблемами.

Я использую jupyter notebook, проблема, с которой я сталкиваюсь, заключается в том, что при вызове nlp.solve (obj_poss) тогда,

1.) Функция solve () не выводит информацию о том, что происходит внутри, как это должно происходить, см.:https://www.coin-or.org/Ipopt/documentation/node36.html

2.) При наличии ошибок в аргументах, отправленных в solve(), не выводится причина возникновения ошибки. Пример попробуйте изменить low_var на low_var = np.array([0.5, 0.5 , 0.5, 0.5]).

Пояснение к коду: Код пытается минимизировать функцию eval_Psum. Под ограничением функции eval_g_func.

 from dolfin import * # Only here to give pyipopt the correct petsc_comm_world
import numpy as np
import pyipopt
import sympy
from pdb import set_trace



def eval_Psum(pos_vec):
    # This function takes in a vector with possitions and calculates the square sum of its
    # coordinates. the vector is aranged as [x1, y1, x2, y2,...]
    sq_sum = 0
    for i in range(len(pos_vec)):
        sq_sum = sq_sum   pos_vec[i]**2
    return sq_sum

def eval_dPsum(pos_vec):
    return np.array([2*pos_vec[0], 2*pos_vec[1], 2*pos_vec[2], 2*pos_vec[3]])

def eval_g_func(pos):
    num_obj = 2
    x = list(sympy.symbols("x0:%d" % num_obj))
    y = list(sympy.symbols("y0:%d" % num_obj))


    g, z = [], []
    for i in range(num_obj):
        z.append(x[i])
        z.append(y[i])
        for j in range(i 1, num_obj):
            xi, yi = x[i], y[i]
            xj, yj = x[j], y[j]
            int_radius = 10**-5
            # FIXME: Max_int_radius should change with each cable radius
            g.append(int_radius**2 - (xi - xj)**2 - (yi - yj)**2)


    g = sympy.Matrix(g)
    z = sympy.Matrix(z)
    print("symbolic g ", g)

    eval_g = g.subs([(z[i], pos[i]) for i in range(2*num_obj)]) 
    return np.array(eval_g).astype(float)



def eval_jac_g_func(pos, flag):
    if flag:
        nvar = len(pos)
        nobj = int(nvar/2)
        ncon = int(nobj   nobj*(nobj-1)/2)
        rows = []
        for i in range(ncon):
            rows  = [i]*nvar
        cols = list(range(nvar))*ncon
        return (np.array(rows), np.array(cols))
    else:

        num_obj = 2
        x = list(sympy.symbols("x0:%d" % num_obj))
        y = list(sympy.symbols("y0:%d" % num_obj))

        g, z = [], []
        for i in range(num_obj):
            z.append(x[i])
            z.append(y[i])
            for j in range(i 1,num_obj):
                xi, yi = x[i], y[i]
                xj, yj = x[j], y[j]
                int_radius = 10**-5
                # FIXME: Max_int_radius should change with each cable radius
                g.append(int_radius**2 - (xi - xj)**2 - (yi - yj)**2)


        g = sympy.Matrix(g)
        z = sympy.Matrix(z)


        jac_g = g.jacobian(z) # Create jacobian of the constraints
        print("symbolic jac_g ", jac_g)

        eval_jac_g = jac_g.subs([(z[i], pos[i]) for i in range(2*num_obj)]) 
        eval_jac_g = np.array(eval_jac_g).astype(float)
        return eval_jac_g





obj_poss = np.array([0.9, 0.9, 0.8, 0.8])
print(eval_jac_g_func(obj_poss, False ) )


nvar = 4 # number of variables x1 x2 x3 x4
low_var = np.array([0, 0,0, 0]) # lower constraint on variables
#low_var = np.array([0.5, 0.5, 0.5, 0.5]) # Uncheck this to see what happends
up_var = np.array([1, 1, 1, 1]) # upper constraint on variables
ncon = 1 # Later we add constraints that turbines cannot overlap # number of constraining equations = 1 namely g 
gl =  - np.inf*np.ones(ncon, dtype=float) # lower constraint value on the constraint function g
gu = np.zeros(ncon, dtype=float) # upper constraint value on the constraint function g
num_non_zero_j = 0 # Number of non zero terms in jacobian
num_non_zero_H = 0 # Number of non zero terms in Hessian


nlp = pyipopt.create(nvar,      # Number of controls
                                  low_var,  # Lower bounds for Control
                                  up_var,   # Upper bounds for Control
                                  ncon,      # Number of constraints
                                  gl,  # Lower bounds for contraints
                                  gu,   # Upper bounds for contraints
                                  nvar*ncon, # Number of nonzeros in cons. Jac
                                  0,         # Number of nonzeros in cons. Hes
                                  lambda x: eval_Psum(x),  # Objective eval
                                  lambda x: eval_dPsum(x), # Obj. grad eval
                                  eval_g_func,    # Constraint evaluation
                                  eval_jac_g_func # Constraint Jacobian evaluation

                         )

obj_poss = np.array([0.9, 0.9, 0.8, 0.8])
nlp.int_option("print_level",12)


solution = nlp.solve(obj_poss)
  

Попробовал быструю репликацию, но не могу найти правильный pyipopt для моей среды. Если они используют logger для вывода сообщений, вам может потребоваться инициализировать logger на уровне «DEBUG», чтобы все было записано в консоль Jupyter, например:

 import logging
logger = logging.getLogger()
logger.setLevel(logging.DEBUG)