我将制作一个脚本来优化潮汐涡轮机在溪流中的位置。为此，我将找到使涡轮机阻力最大化的位置。
我打算使用 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)