How does `KSComp` behave when using the options `lower_flag`, `ref`, `ref0` and `units` and `upper` at the same time?

Reading the Adding Constraints page, I have understood that the scaling of the a constraint value happens in the following order:

the constraint value is converted from the connected output’s units to the desired units;
ref0 and ref scale the constraint value to 0 and 1, respectively, and they have to be defined in the desired units.

My understanding is that also lower and upper have to defined in the desired units and that they are scaled according to the values of ref0 and ref. Is this correct?

Now, I'm finding a little bit trickier to understand what happens with KSComp. Say I have the following situation:

the connected output feeds an array of stiffnesses defined in N/mm;
the KSComp component is defined with the following options: lower_flag=True, units=N/m and upper=1, meaning that the stiffnesses will be defined in N/m and will be constrained to be larger than 1.

My understanding is that the array of stiffnesses will be converted to N/m and then the component will transform the input array and calculate the KS function such that the constraint is satisfied if it is less than zero.

How should I use ref0 and ref in this case? In theory they should represent the physical values in the desired units that map to 0 and 1, respectively. However, looking at the source code of KSComp it seems that the scaling is applied to the computed KS function. In this case, if I only use the ref option and I scale a feasible value to 1, this will never satisfy the constraint, because with the option add_constraint=True the constraint is set to be satisfied when the KS function si smaller than 0 (0 scales to 0). Am I getting any of this wrong?

Solution

When adding a constraint the values of lower, upper, scaler, adder, ref and ref0 should all be applied in the units of the constraint.

The KS function behind KS comp can be a bit difficult to scale because it produces an output that is somewhat nonintuitive.

I massaged the add_constraint example from the docs to make upper easily changed, and had the example just run the model rather than optimize it.

It prints out the quantity KS being produced (and constrained) by the KSComp. This quantity is positive when violated, so 0 is probably a good value for ref0 (which is the default. Playing around with values for upper indicates that the nominal value of KS is about the same order of magnitude of the distance from the constraint. So if your initial guess has a minimum distance of 3000 units from upper, the out KS will be on the same order of magnitude, 3000.

If we assume to scale this so that the nominal constraint values have a magnitude of 1 (not always the best assumption but a reasonable one), then an appropriate value for ref in this case is 3000.

(I think it's always best to let ref be positive).

import numpy as np
import matplotlib.pyplot as plt
import openmdao.api as om

n = 50
upper = -40
prob = om.Problem()
model = prob.model

prob.driver = om.ScipyOptimizeDriver()

model.add_subsystem('ivc', om.IndepVarComp('x', val=np.ones(n,), units='N/mm'))

model.add_subsystem('comp', om.ExecComp('y = 3.0*x**2 + k',
                                        x={'val': np.zeros((n, )), 'units': 'N/m'},
                                        y={'val': np.zeros((n, )), 'units': 'N/m'},
                                        k=0.0), promotes_inputs=['x', 'k'])

model.add_subsystem('ks', om.KSComp(width=n, lower_flag=True, upper=upper,
                                    ref=40, add_constraint=True))

model.add_design_var('k', lower=-100, upper=100)
model.add_objective('k', scaler=1)

model.connect('comp.y', 'ks.g')

prob.setup()
prob.set_val('x', np.linspace(-np.pi/2, np.pi/2, n))
prob.set_val('k', 5.)

prob.run_driver()

fig, ax = plt.subplots()

x = prob.get_val('x')
y = prob.get_val('comp.y')
ks = prob.get_val('ks.KS')

print(f'KS = {ks}')

ax.plot(x, y, 'r.')
ax.plot(x, upper*np.ones_like(x), 'k--')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.grid(True)
ax.text(-0.25, 0, f"k = {prob.get_val('k')[0]:6.3f}")

plt.show()