basf / doe Goto Github PK

David Hajnal recently asked for adding support for nonlinear constraints. This issue addresses this request.

NchooseK constraint fails when in addition a mixture constraint (sum(x_i) = 1.0) is included. The warning Sampling of points fulfilling this problem's constraints is not implemented appears. The resulting design is a trivial collection of runs with a single input set to 1 with all the other to 0. In the example code below, the goal is to have a simple linear design for a four-ingredient mixture (x1+x2+x3+x4+x5=1) with a maximum of three ingredients:

num_inputs = 4
inputs = [f"x{i+1}" for i in range(num_inputs)]
print('Inputs\t   :', inputs)

# mixture constraint
constr1 = opti.LinearEquality(names=[x for x in inputs], lhs=[1 for _ in inputs], rhs=1)
print('Constraints:', constr1)
# n-choose-k
num_nCk_active = 2
constr2 = opti.NChooseK([x for x in inputs], max_active=num_nCk_active)

problem = opti.Problem(
    inputs=[opti.Continuous(x, domain=[0., 1.]) for x in inputs],
    outputs=[opti.Continuous("y")],
    constraints=[constr2, constr1],
)

doe = doe.find_local_max_ipopt(
        problem=problem,
        model_type="linear"
)

A nice-to-have method could return a D-optimality score, such as the log-determinant of the information matrix log(det(X^T X)). This would allow easy comparison between designs.

There are currently no tests for the ProblemContext class.

DoE depends on opti and at this point essentially provides a single user-facing function.
I propose to merge doe into opti. Since ipopt is not a nicely pip-installable depency, this should be done as an optional submodule.

The advantage would be that generic functionality such as the relaxation problem context and the corner point sampling would profit opti as well, and we have less hassle developing a single package.

If no jacobian is provided for a constraint, IPOPT will approximate it using finite differences. For linear constraints the jacobian is generated internally by doe during the processing of the provided opti.Problem. For arbitrary nonlinear constraints this cannot be done easily. However, in high dimensions finite the difference approximation is too time consuming.
Therefore, there should be an option for the user to provide a jacobian for nonlinear constraints.

To illustrate the problem: Here we have 10 variables with 1 nonlinear constraint and a linear model.

problem = opti.Problem(
    inputs=opti.Parameters([opti.Continuous(f"x{i+1}", [0, 1]) for i in range(10)]),
    outputs=[opti.Continuous("y")],
    constraints=[opti.NonlinearEquality("(x1**2 + x2**2 + x3**2 + x4**2 + x5**2 + x6**2 + x7**2 + x8**2 + x9**2)**0.5 - x10")],
)

result = find_local_max_ipopt(problem, "linear", ipopt_options={"disp":5, "maxiter": 100})
result.round(3)

And now the same with a linear constraint (mixture constraint)

problem = opti.Problem(
    inputs=opti.Parameters([opti.Continuous(f"x{i+1}", [0, 1]) for i in range(10)]),
    outputs=[opti.Continuous("y")],
    constraints=[opti.LinearEquality(names=[f"x{i+1}" for i in range(10)], rhs=1)],
)

result = find_local_max_ipopt(problem, "linear", ipopt_options={"disp":5, "maxiter": 100})
result.round(3)

When you run both examples you will see that even in case of a moderate problem size of 10 variables with only one nonlinear constraint we observe that the optimization slows down significantly (and maybe even get worse results due to inaccuracies of the finite difference approximation).

When a design does not satisfy a constraint it is returned as if everything went well, it would be nice to warn the user about this.
Example: a simple three-input linear design with a formulation constraint (x1+x2+x3=1). If all inputs have the domain [0.5,1], then the formulation constraint cannot be satisfied

import pandas as pd
import opti
import doe

inputs = ['x1', 'x2', 'x3']
# formulation constraint
constr1 = opti.LinearEquality(names=inputs, rhs=1)
problem = opti.Problem(
                      inputs=[opti.Continuous(x, domain=[0.5, 1]) for x in inputs],
                      outputs=[opti.Continuous("y")],
                      constraints=[constr1],
)
df = doe.find_local_max_ipopt(problem=problem,
                              model_type="linear")

constr1.satisfied(df).values
# array([False, False, False, False, False, False, False])

The last output shows that the formulation constraint is not satisfied; nonetheless, the (incorrect) design is generated and returned with no warnings to the user.

We currently temporarily increase python's builtin recursion limit as a workaround for matthewwardrop/formulaic#78.
With formulaic 0.5 we can remove that workaround.

The method find_local_max_ipopt could include the possibility of including center points and/or replicates in the design.