I am trying to translate a Matlab code to python, and it includes a multi-objective solver (using gamultiobj in matlab) for which I chose to use the AGEMOEA algorithm in the pymoo-package. However, my problem is that I can't seem to understand how to format the input correctly, and as a result, the pymoo-solver is not able to reshape the output.
My error is this:
File "C:\Users\IDABUK\Anaconda3\lib\site-packages\pymoo\core\problem.py", line 288, in _format_dict
raise Exception(
Exception: ('Problem Error: F can not be set, expected shape (200, 1) but provided (200, 1, 1, 1, 2)', ValueError('cannot reshape array of size 400 into shape (200,1)'))
And here is the code (I tried to minimize it)
Imports and defining variables
import numpy as np
import warnings
from pymoo.core.problem import Problem
from pymoo.algorithms.moo.age import AGEMOEA
from pymoo.optimize import minimize
from pymoo.problems.functional import FunctionalProblem
from pymoo.termination import get_termination
from cmath import nan
import numpy as np
from scipy.optimize import fsolve
rho = 1000
mu = 9 * 10 ** - 4
roughness = 0.2 * 0.001
D_min = 5 * 0.001
D_max = 1000 * 0.001
velocity_boundaries = np.array([0.1, 3])
m_dot_plot = np.linspace(0.1, 1, 10)
D_plot = np.array([15, 20, 30, 40, 50, 65, 80, 100, 125, 150, 200, 250,
300, 350, 400]) / 1000
Here is the code
def colebrook(Re=None, k=None):
if Re <= 0 or k <= 0:
raise Exception('inputs must be all positive values')
else:
if Re < 2300:
f = 64 / Re
message = 1
else:
def fun(f):
return 1 / np.sqrt(f) + 2 * np.log10(k / 3.7 + 2.51 / (Re * np.sqrt(f)))
f, _, flag, msg = fsolve(fun, 0.1, full_output=True) # removed optimset('Display', 'off') #ERROR! FIXME
# discuss solutions
if flag <= 0:
# refuse solution
f = nan
message = - 1
# algorithm did not converge to a solution.
else:
# accept solution
message = 1
return f, message
def nlinconstraint(D=None, m_dot=None, rho=None, velocity_boundaries=None):
v = 4 * m_dot / (np.pi * D ** 2.0 * rho)
# constraint functions
if True and not len(velocity_boundaries) == 0:
if np.isscalar(velocity_boundaries[0]) and not np.isnan(velocity_boundaries[0]): #wrong with () instead of []? XXX
#c[:, 0] = velocity_boundaries[0] - v
c1= np.array(velocity_boundaries[0] - v)
if np.isscalar(velocity_boundaries[1]) and not np.isnan(velocity_boundaries[1]):
#c[:, 1] = v - velocity_boundaries[1]
c2 = np.array(v - velocity_boundaries[1])
#c = np.array([c1, c2])
c = np.concatenate((c1, c2), axis= 0)
ceq = np.array([])
return c, ceq #remove ceq as no equalities as constraints?
def obj_function(D=None, m_dot=None, rho=None, mu=None, roughness=None):
# This function contains the objectives to be minimized.
# The first objective is expressed as the pressure drop per meter of
# pipe, i.e. [Pa/m]. The second one is the inner pipe diameter in [m]
# To run this function, the Colebrook implicit solver is required.
#
# pressure drop [Pa/m]
# delta p 8 * m^2 4 m epsilon
# -------- = ----------------- * colebrook(---------- , --------- )
# L pi^2 * rho * D^5 pi * D * m D
#y = np.array([]) # empty array
v = 4 * m_dot / (np.pi * D ** 2.0 * rho)
Re = (np.multiply(np.multiply(rho, v), D)) / mu
f = np.array([])
f = colebrook(Re, roughness/D)[0] # do not include message?
delta_p_m = np.multiply(np.multiply(np.multiply((1.0/D), f) * 0.5, rho),
v**2)
D = np.array(D)
y = np.column_stack((delta_p_m, D))
return y
def topsis(DM=None, weights=None, min_max=None):
# This function ranks the decision matrix [m x n]
# with n = attributes (columns)
# and m = alternatives, configurations (rows) using the TOPSIS method.
#
# Inputs: - decision matrix [m x n]
# - weights array for the attributes [1 x m]
# - minimize/maximize array [1 x m]
# for the j-th attribute, with j = 1, 2, ... m
# choose +1 if the j-th attribute shall be maximized
# -1 if the j-th attribute shall be minimized
#
# Output: - ranked decision matrix [m x n]
# expressed with the "real" units.
#
# Pay attention, this function works if each element of the decision matrix
# is a scalar, not an array!
m, n = DM.shape # fix
Z = np.multiply(nan, np.ones((m, n)))
Zw = Z
a_plus = np.multiply(nan, np.ones((1, n)))
a_minus = np.multiply(nan, np.ones((1, n)))
for j in np.arange(1, n+1).reshape(-1): # fix
Z[:, j] = DM[:, j]/(np.sum(DM[:, j]**2, 1-1))**0.5
Zw[:, j] = np.multiply(Z[:, j], weights(j))
if min_max(j) == -1: # is this a typo? Array index or function call?
a_plus[j] = np.amin[Zw[:, j]]
a_minus[j] = np.amax[Zw[:, j]]
else:
if min_max(j) == 1: # same as above
a_plus[j] = np.amax[Zw[:, j]]
a_minus[j] = np.amin[Zw[:, j]]
# 2-1 for indexing purposes
Si_minus = np.sum((Zw - a_minus)**2, 2-1)**0.5
Si_plus = np.sum((Zw - a_plus) ** 2, 2-1)**0.5
c = Si_minus/(Si_plus + Si_minus)
new_DM = np.array([DM, c])
# ranked_DM = sortrows(new_DM,n + 1,'descend') #fix - matlab code
ranked_DM = new_DM.sort(reverse=True) # fix - include n+1 from matlab code -use argsort instead (FIXFIX) # noqa: E501
output_DM = ranked_DM[:, np.arange(1, n+1)]
return output_DM
def opt_pr_dia(m_dot=None, rho=None, mu=None, roughness=None,
D_min=None, D_max=None, velocity_boundaries=None):
# This function solves the optimization problem of the pipe optimal
# diameter with minimal pressure drop.
# Inputs: - mass flow rate [kg/s] (scalar)
# - fluid density [kg/m^3] (scalar)
# - pipe absolute roughness [m] (scalar)
# - lower boundary for the search of the pipe D [m] (scalar)
# - upper boundary for the search of the pipe D [m] (scalar)
# - velocity boundary [v_min, v_max] (1 x 2 array)
#
# Outputs: - T table containing all non-dominated solutions
# - exit flag (>0 if successful, <=0 if algorithm did not
# converge)
f_array = lambda D = None: [obj_function(D,m_dot,rho,mu,roughness)]
constr_ieq = [lambda D = None: nlinconstraint(D, m_dot, rho, velocity_boundaries)]
n_var = 1 #number of variables
problem = FunctionalProblem(n_var, f_array, constr_ieq=constr_ieq, xl=D_min, xu=D_max ) #BUG: gets shape of problem (200, 1, 1, 2, 1) instead of (200, 1) - ValueError cannot reshape array of size 400 into shape (200,1) - size is twice what it should be
#problem = Problem(n_var, f_array, n_ieq_constr=constr_ieq, xl=D_min, xu=D_max)
algorithm = AGEMOEA(pop_size= 200) #choice of algorithm with chosen population size
termination = get_termination("n_gen", 500) # terminates after 500 generations - could also terminate due to time: ("time", "00:05:00") - 5 minutes
res = minimize(problem, algorithm, termination, seed=1, verbose=True)
f_values = algorithm.result()
exitflag = 500 - res.algorithm.n_gen #number of iterations left
if exitflag < 0:
# refuse solution
raise Exception('optimization failed')
else:
if exitflag == 0:
# max number of generations reached
# FIXME: Implement using Python logging - done(?)
warnings.warn('max number of generations reached. ' +
'Try to relax variable boundaries or constraints')
else:
# accept solution
print('optimization successful')
w = np.array([0.5, 0.5])
# rank the objectives according to the TOPSIS method
obj = topsis(f_values, w, np.array([-1, -1]))
# build the output table containing the designs in the pareto front
d = obj[:, 2]
delta_p_1m = obj[:, 1]
v = 4 * m_dot/(np.pi * d**2.0 * rho)
Re = np.multiply(np.multiply(rho, v), d)/mu
delta_p_mmWC = delta_p_1m/9.81
T = np.array(np.multiply(m_dot, np.ones((delta_p_1m.shape, # does this work??
delta_p_1m.shape))),
np.round(delta_p_1m), np.round(delta_p_mmWC),
np.round(d * 1000, 2), np.round(v, 2), np.round(Re))
#FIXME??
T.Properties.VariableNames = np.array(['m_dot', 'dp [Pa/m]',
'dp [mm H2O/m]', 'd [mm]',
'v [m/s]', 'Re'])
return T, exitflag
if __name__ == '__main__':
# test the function
m_dot = 0.5
rho = 1000
mu = 0.001
roughness = 0.00005
D_min = 0.01
D_max = 0.05
velocity_boundaries = np.array([0.5, 2])
T, exitflag = opt_pr_dia(m_dot, rho, mu, roughness, D_min, D_max, velocity_boundaries)
Here is the original code in matlab: https://github.com/giuliots95/Optimize_pipe_diameter_and_pressure_drop
