diff --git a/src/gemseo/problems/optimization/hock_schittkowski_71.py b/src/gemseo/problems/optimization/hock_schittkowski_71.py new file mode 100644 index 0000000000000000000000000000000000000000..9d9c8c0d255c520bf999ab952cd40bc8fd736d58 --- /dev/null +++ b/src/gemseo/problems/optimization/hock_schittkowski_71.py @@ -0,0 +1,211 @@ +# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com +# +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License version 3 as published by the Free Software Foundation. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program; if not, write to the Free Software Foundation, +# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +r"""Hock & Schittkowski problem 71. + +This module implements the Hock & Schittkowski non-linear programming problem 71. + +See: Willi Hock and Klaus Schittkowski. (1981) Test Examples for +Nonlinear Programming Codes. Lecture Notes in Economics and Mathematical +Systems Vol. 187, Springer-Verlag. +Based on MATLAB code by Peter Carbonetto. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +from numpy import array + +from gemseo.algos.design_space import DesignSpace +from gemseo.algos.optimization_problem import OptimizationProblem +from gemseo.core.mdo_functions.mdo_function import MDOFunction + +if TYPE_CHECKING: + from gemseo.typing import NumberArray + + +class HockSchittkowski71(OptimizationProblem): + """Hock and Schittkowski problem 71.""" + + def __init__( + self, + initial_guess: NumberArray = (1.0, 5.0, 5.0, 1.0), + ): + """Initialize the Hock Schittkowski 71 problem. + + Args: + initial_guess: The initial guess for the optimal solution. + """ + design_space = DesignSpace() + design_space.add_variable( + "x", + size=4, + lower_bound=[1.0, 1.0, 1.0, 1.0], + upper_bound=[5.0, 5.0, 5.0, 5.0], + ) + + design_space.set_current_value(array(initial_guess)) + + super().__init__(design_space) + + self.objective = MDOFunction( + self.compute_objective, + name="hock_schittkoski_71", + f_type=MDOFunction.FunctionType.OBJ, + jac=self.compute_objective_jacobian, + expr="x_1 * x_4 * (x_1 + x_2 + x_3) + x_3", + input_names=["x"], + dim=1, + ) + + equality_constraint = MDOFunction( + self.compute_equality_constraint, + name="equality_constraint", + f_type=MDOFunction.ConstraintType.EQ, + jac=self.compute_equality_constraint_jacobian, + expr="(x_1**2 + x_2**2 + x_3**2 + x_4**2) - 40", + input_names=["x"], + dim=1, + ) + self.add_constraint( + equality_constraint, constraint_type=MDOFunction.ConstraintType.EQ + ) + + inequality_constraint = MDOFunction( + self.compute_inequality_constraint, + name="inequality_constraint", + f_type=MDOFunction.ConstraintType.INEQ, + jac=self.compute_inequality_constraint_jacobian, + expr="25 - (x_1 * x_2 * x_3 * x_4)", + input_names=["x"], + dim=1, + ) + self.add_constraint( + inequality_constraint, constraint_type=MDOFunction.ConstraintType.INEQ + ) + + @staticmethod + def compute_objective(design_variables: NumberArray) -> NumberArray: + """Compute the objectives of the Hock and Schittkowski 71 function. + + Args: + design_variables: The design variables vector. + + Returns: + The objective function value. + """ + return ( + design_variables[0] * design_variables[3] * (sum(design_variables[0:3])) + + design_variables[2] + ) + + @staticmethod + def compute_objective_jacobian(design_variables: NumberArray) -> NumberArray: + """Compute the Jacobian of objective function. + + Args: + design_variables: The design variables vector. + + Returns: + The gradient of the objective functions wrt the design variables. + """ + x_1, x_2, x_3, x_4 = design_variables + + return array([ + x_1 * x_4 + x_4 * (x_1 + x_2 + x_3), + x_1 * x_4, + x_1 * x_4 + 1.0, + x_1 * (x_1 + x_2 + x_3), + ]) + + @staticmethod + def compute_equality_constraint(design_variables: NumberArray) -> NumberArray: + """Compute the equality constraint function. + + Args: + design_variables: The design variables vector. + + Returns: + The equality constraint's value. + """ + x_1, x_2, x_3, x_4 = design_variables + return array([x_1**2 + x_2**2 + x_3**2 + x_4**2 - 40]) + + @staticmethod + def compute_inequality_constraint(design_variables: NumberArray) -> NumberArray: + """Compute the inequality constraint function. + + Args: + design_variables: The design variables vector. + + Returns: + The inequality constraint's value. + """ + x_1, x_2, x_3, x_4 = design_variables + return array([25 - x_1 * x_2 * x_3 * x_4]) + + @staticmethod + def compute_equality_constraint_jacobian( + design_variables: NumberArray, + ) -> NumberArray: + """Compute the equality constraint's Jacobian. + + Args: + design_variables: The design variables vector. + + Returns: + The Jacobian of the equality constraint function wrt the design variables. + """ + x_1, x_2, x_3, x_4 = design_variables + + return array([ + 2 * x_1, + 2 * x_2, + 2 * x_3, + 2 * x_4, + ]) + + @staticmethod + def compute_inequality_constraint_jacobian( + design_variables: NumberArray, + ) -> NumberArray: + """Compute the inequality constraint's Jacobian. + + Args: + design_variables: The design variables vector. + + Returns: + The Jacobian of the inequality constraint function wrt the design variables. + """ + x_1, x_2, x_3, x_4 = design_variables + + return array([ + -x_2 * x_3 * x_4, + -x_1 * x_3 * x_4, + -x_1 * x_2 * x_4, + -x_1 * x_2 * x_3, + ]) + + @staticmethod + def get_solution() -> tuple[NumberArray, NumberArray]: + """Return the analytical solution of the problem. + + Returns: + The theoretical optimum of the problem. + """ + x_opt = array([0.99999999, 4.74299964, 3.82114998, 1.37940829]) + f_opt = x_opt[0] * x_opt[3] * sum(x_opt[0:3]) + x_opt[2] + return x_opt, f_opt diff --git a/tests/problems/optimization/test_analytical.py b/tests/problems/optimization/test_analytical.py index 834c9f18fa7d9fa5a8db90b7bc4a8ba00df0c154..d0cd52bdc69a5ea8d7993a8f8f13ddf10be2381a 100644 --- a/tests/problems/optimization/test_analytical.py +++ b/tests/problems/optimization/test_analytical.py @@ -27,6 +27,7 @@ from numpy import zeros from gemseo.algos.base_driver_library import MaxIterReachedException from gemseo.algos.opt.factory import OptimizationLibraryFactory +from gemseo.problems.optimization.hock_schittkowski_71 import HockSchittkowski71 from gemseo.problems.optimization.power_2 import Power2 from gemseo.problems.optimization.rastrigin import Rastrigin from gemseo.problems.optimization.rosen_mf import RosenMF @@ -83,3 +84,9 @@ def test_rosen_mf() -> None: step=1e-8, threshold=1e-4, ) + + +def test_hs071() -> None: + """""" + problem = HockSchittkowski71() + run_and_test_problem(problem) diff --git a/tests/problems/optimization/test_hock_schittkowski_71.py b/tests/problems/optimization/test_hock_schittkowski_71.py new file mode 100644 index 0000000000000000000000000000000000000000..ad49658b514b19335eecfa54c205e6a330c1c243 --- /dev/null +++ b/tests/problems/optimization/test_hock_schittkowski_71.py @@ -0,0 +1,114 @@ +# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com +# +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License version 3 as published by the Free Software Foundation. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program; if not, write to the Free Software Foundation, +# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +"""Tests for the Hock and Schittkowski problem 71.""" + +from __future__ import annotations + +import pytest +from numpy import array + +from gemseo.problems.optimization.hock_schittkowski_71 import HockSchittkowski71 + + +@pytest.fixture +def hock_schittkowski_71() -> HockSchittkowski71: + """Create a :class:`.HockSchittkowski71` optimization problem. + + Returns: + A HockSchittkowski71 instance. + """ + return HockSchittkowski71(initial_guess=array([1.0, 5.0, 5.0, 1.0])) + + +def test_obj_jacobian(hock_schittkowski_71): + """Test the Jacobian of the Hock and Schittkowski 71 objective function. + + Args: + hock_schittkowski_71: Fixture returning a HockSchittkowski71 + `OptimizationProblem`. + """ + x_dv = array([1.0, 1.0, 1.0, 1.0]) + hock_schittkowski_71.objective.check_grad(x_dv, error_max=1e-6) + + +def test_equality_constraint_jacobian(hock_schittkowski_71): + """Test the equality constraint function's Jacobian of Hock and Schittkowski 71. + + Args: + hock_schittkowski_71: Fixture returning a HockSchittkowski71 + `OptimizationProblem`. + """ + x_dv = array([1.0, 1.0, 1.0, 1.0]) + hock_schittkowski_71.constraints[0].check_grad(x_dv, error_max=1e-6) + + +def test_inequality_constraint_jacobian(hock_schittkowski_71): + """Test the inequality constraint function's Jacobian of Hock and Schittkowski 71. + + Args: + hock_schittkowski_71: Fixture returning a HockSchittkowski71 + `OptimizationProblem`. + """ + x_dv = array([1.0, 1.0, 1.0, 1.0]) + hock_schittkowski_71.constraints[1].check_grad(x_dv, error_max=1e-6) + + +def test_compute_objective(hock_schittkowski_71): + """Test the objective function of the Hock and Schittkowski 71 problem. + + Args: + hock_schittkowski_71: Fixture returning a HockSchittkowski71 + `OptimizationProblem`. + """ + x_dv = array([1.0, 1.0, 1.0, 1.0]) + objective = hock_schittkowski_71.compute_objective(x_dv) + + assert objective == 4.0 + + +def test_compute_equality_constraint(hock_schittkowski_71): + """Test the equality constraint function of the Hock and Schittkowski 71 problem. + + Args: + hock_schittkowski_71: Fixture returning a HockSchittkowski71 + `OptimizationProblem`. + """ + x_dv = array([1.0, 1.0, 1.0, 1.0]) + equality_constraint = hock_schittkowski_71.compute_equality_constraint(x_dv) + assert equality_constraint == -36.0 + + +def test_compute_inequality_constraint(hock_schittkowski_71): + """Test the inequality constraint function of the Hock and Schittkowski 71 problem. + + Args: + hock_schittkowski_71: Fixture returning a HockSchittkowski71 + `OptimizationProblem`. + """ + x_dv = array([1.0, 1.0, 1.0, 1.0]) + inequality_constraint = hock_schittkowski_71.compute_inequality_constraint(x_dv) + assert inequality_constraint == 24.0 + + +def test_solution(hock_schittkowski_71) -> None: + """Check the objective value at the solution. + + Args: + hock_schittkowski_71: Fixture returning a HockSchittkowski71 + `OptimizationProblem`. + """ + x_opt, f_opt = hock_schittkowski_71.get_solution() + assert hock_schittkowski_71.objective.evaluate(x_opt) == pytest.approx(f_opt)