Min Max component, objective function

user_16421467 0
Viewed 64

I would like to perform some optimizations by minimizing the maximum of a specific path variable within Dymos. or the maximum of the absolute of such a variable.

In linear programming methods, this can be done by introducing slack variables.

Do you know if this has been attempted before with Dymos, or if there was a reason not to include it?

I understand gradient based methods are not entirely suitable for these problems, though I think some "functions" can be introduced to mitigate this.

For example, The space shuttle reentry problem from [Betts][1] used as a [test example][2] in dymos, the original source contains an example where the maximum heat flux is minimized. Such functionality could be implemented with the "loc" argument as:

phase.add_objective('q_c', loc='max')

[1]: J. Betts. Practical Methods for Optimal Control and Estimation Using Nonlinear Programming. Society for Industrial and Applied Mathematics, second edition, 2010. URL: https://epubs.siam.org/doi/abs/10.1137/1.9780898718577, arXiv:https://epubs.siam.org/doi/pdf/10.1137/1.9780898718577, doi:10.1137/1.9780898718577. [2]: https://openmdao.github.io/dymos/examples/reentry/reentry.html

openmdao
2 Answers

This has been done with pseudospectral methods before. Dymos currently doesn't have any direct way of implementing this, for a few reasons:

  1. As you said, doing this naively can introduce discontinuous gradients that confuse the optimizer. When the node at which the maximum occurs switches, this tends to cause a sharp edge discontinuity in the gradient.

  2. Since the pseudospectral methods are discrete, you cannot guarantee that the maximum will occur at a node. It's often fine to assume it does, but sometimes your requirements might demand more precision.

There are two possible ways to get around this.

  1. The KSComp in OpenMDAO can be used as a "differentiable maximum". Add one after the trajectory, feed it the timeseries data for the output of interest, and set it up such that it returns a smooth approximation to the maximum. The KS function is a bit conservative, so it won't pick out the precise maximum, but depending on the value of the rho option it can be tuned to get pretty close.

  2. When a more precise value of a maximum is needed, it's pretty common to set up a trajectory such that a phase ends when the maximum or minimum is reached.

If the variable whose maximum is being sought is a state, this can be done by adding a boundary constraint on the rate source for that state. This ensures that the maximum occurs at the first or last node in the phase (depending on if its an initial or final boundary constraint). That lets you more accurately capture its value.

If the variable being sought is not a state, its possible to use the polynomials that are used for fitting states and controls in a phase to interpolate the variable of interest. By then taking the time derivative of that polynomial we can get a reasonably good approximation for its rate. The master branch of dymos has a method add_timeseries_rate_output that does this. And soon, within a few weeks hopefully, we'll add add_boundary_rate_constraint so that these interpolated rates can be easily used as boundary constraints.

In the meantime, you should be able to achieve this by adding the timeseries rate output and then manually applying the OpenMDAO method 'add_constraint' to the resulting timeseries output, using either indices=[0] or indices=[-1] to treat it as an initial or final constraint.

This is a common enough request that we'll add some documentation on how to achieve this behavior using both the KSComp approach and the boundary constraint approach.

user_7545360

Personally I'm not as much of a fan of KSComp because I've had trouble getting problems getting those types of objectives to converge in the past. I've used the slack variable and that has worked well. In the following example, we take a guess at the Rotor power in static analysis, and then we run a trajectory and get the actual rotor power during the mission. The objective was to minimize aircraft weight, so if you have a large amount of power in statics, that costs more weight. The constraint shown below prevents us from decreasing our updated guess of rotor power in statics below the maximum power required during the trajectory.

p.model.add_subsystem(
       'static_power_check', 
       om.ExecComp('Power_check = Power_ODE - Power_statics',
                   Power_check = {'value':np.ones(nn_timeseries_main_tx), 'units':'kW'},
                   Power_ODE = {'value':np.ones(nn_timeseries_main_tx), 'units':'kW'},
                   Power_statics = {'value':0.0, 'units':'kW'}),
       promotes_inputs=[
               ('Power_ODE','hop0.main_phase.timeseries.Power_R'), ('Power_statics','Power_{rotor,slack}')],
       promotes_outputs=['Power_check'])
p.model.add_constraint('Power_check', upper=0, ref=1)

The constraint on the slack variable effectively helped us ensure that our slack rotor power matched the maximum rotor power during the mission. This allowed us to get the right sizes for the rotor parts (i.e. motors).

user_198963780
Related