Pymoo - 无法重塑数组，求解器输入错误（matlab 到 python 翻译）

问：

我正在尝试将 Matlab 代码转换为 python，它包括一个多目标求解器（在 matlab 中使用 gamultiobj），为此我选择在 pymoo-package 中使用 AGEMOEA 算法。 但是，我的问题是我似乎无法理解如何正确格式化输入，因此，pymoo-solver 无法重塑输出。

我的错误是这样的：

     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)'))

这是代码（我试图最小化它）

导入和定义变量

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

这是代码

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)

以下是matlab中的原始代码： https://github.com/giuliots95/Optimize_pipe_diameter_and_pressure_drop