Tuesday, June 4, 2019

Logarithmic mean temperature difference requires yet another cone?

The logarithmic mean temperature difference
$$LMTD(x,y) = \frac{x-y}{\ln(x/y)}$$

can be used to model heat exchange, e.g., when optimizing the reuse of excess process heat onsite at industrial facilities. In the heat exchanger network synthesis problem due to Yee and Grossmann (link goes to a recent paper by Mistry and Misener), we find that its reciprocal raised to some power, that is
$$RecLMTD^\beta(x,y) = \left(\frac{\ln(x/y)}{x-y}\right)^\beta,$$

can be extracted as a separately contributing term in the objective function. Capitalizing on the convexity of this term on $(x,y) \in \mathbb{R}_+^2$, for all considered $\beta \geq 0$, this leads to better performance when solving the otherwise nonconvex problem as argued in the paper.

A challenge to find the conic reformulation of this function was posed under the Oberwolfach Workshop on Mixed-Integer Nonlinear Optimization (2019) and we accepted. Of course, this is trivial if no restrictions are put on the set of cones as one may just define
$$\mathcal{K} = \mathrm{cl}\{(t,s,x,y) \in \mathbb{R}_{++}^4 : t \geq s \cdot RecLMTD^\beta(x/s, y/s) \}$$
and call it a day. This cone is nonempty, closed and convex and hence obeys $K = (K^{*})^{*}$ as well as all the usual properties of conic duality. Computationally, however, the cone is not particularly desirable and we can do better with a bit of reformulation:
$$\begin{array}{lll} t \geq \left(\frac{\ln(x/y)}{x-y}\right)^\beta,\\t \geq \left(\frac{\ln(u/y + 1)}{u}\right)^\beta,& u = x-y,\\ y \geq \frac{u}{\exp(ut^{1/\beta})-1},& u = x-y,\\ y \geq \frac{u}{\exp(u/s)-1},& u = x-y, & s \geq t^{-1/\beta},\\ \end{array}$$
where I substitute in the first step, rewrite assuming either $u > 0$ or $u <0$ (both leads to the same) in the second, and extract a power cone representable subexpression in the third. This means that the representation problem of $RecLMTD^\beta$ have been reduced to the representation problem of
$$\mathcal{K} = \mathrm{cl}\left\{(t,s,x) \in \mathbb{R}_{+}^2 \times \mathbb{R}_{++} : t \geq \frac{x}{\exp(x/s)-1} \right\},$$
which, just like the quadratic, power and exponential cones, is defined as the epigraph of the perspective of a univarite convex function; in this case $\frac{x}{\exp(x)-1}$. Whether this cone can be written in terms of the others, or has potential for computationally efficient implementations itself, remains open. We invite anyone interested in barrier functions and interior-point algorithms to take a crack at it.

What is possible right now?
The logarithmic mean is bounded from above and below by, respectively, the arithmetic and the geometric mean, both of which are representable in MOSEK. There is also Chen's approximation which corresponds to a power cone. Finally it can be observed that $\log(\exp(-x)+1)$ is a fairly good underestimator of $\frac{x}{\exp(x)-1}$. This means that the set described by the inequality from above,
$$y \geq \frac{u}{\exp(u/s)-1},$$

can be outer approximated by
$$y \geq s \log(\exp(-u/s)+1),$$

which is representable in MOSEK as the homogenization of the softplus function. Beware that the usefulness of any of these alternatives is unknown to the author.

Monday, May 20, 2019

Beyond 3/2 - power cones from user perspective

Here is feedback on power cones in Mosek 9 we received (published with permission of the user, anonymous for the purpose of this post).

Back to the power cones: we were excited about this feature because it allows us to implement market impact constraints with exponents that are arbitrary greater than one.

We benchmarked a simple portfolio optimization problem. The goal was to maximize an expected return on a 10 000 instrument portfolio under a market impact constraint. We chose this particular problem because it contains ten thousand power cones in order to implement the exponent part of the market impact constraint.

With Mosek 8.1 we were limited to powers 1, 3/2, 5/3 and 2 because other powers would have taken a gigantic effort to implement. With Mosek 9 all the powers above one are now easily available. As can be seen in the following timings (in seconds), in our case, the power cones are roughly as fast as the implementation using multiple rotated quadratic cones. This is great news as we were expecting a 10x degradation. Some oddities though: power 5/3 seems to be slower in Mosek 9 than in Mosek 8.1 and power 2 seemed faster in power cones than in rotated cones.

VersionMosek 8.1Mosek 9.0Mosek 9.0

ImplementationRotated ConesRotated ConesPower Cones

In conclusion, we are really happy with Mosek 9. The power cones are going to replace all the rotated cone shenanigans we previously had, they are as fast and a lot more powerful/simpler to use.

Tuesday, May 7, 2019

Mosek in Pyomo and Wasserstein distances

Pyomo is an open-source package for modeling optimization problems. You can read about its capabilities here.

The latest release of Pyomo supports MOSEK as a solver for linear, quadratic and quadratically constrained convex problems and their mixed-integer versions. This is thanks to a plugin prepared by our student employee Harun Kivril and thanks to the efforts of the Pyomo development team.

On this occasion we made two notebooks showcasing the use of our own Fusion API, of CVXPY and of Pyomo to model the same problem, namely finding the Wasserstein barycenter of discrete probability distributions, or, in our example application, averaging a set of images of the same handwritten digit.

The first notebook implements a linear Wasserstein barycenter problem in Fusion, CVXPY and Pyomo:

The second one implements an entropy-regularized version of the problem using the exponential cone capability available through Fusion and CVXPY:

Thursday, April 25, 2019

Response codes and statuses

After optimizing a problem with MOSEK it is time to retrieve the results. This post explains the philosophy behind the responses and statuses MOSEK returns.

The first thing to know is if the optimization completed successfully. This is indicated by the response code. Anything other than "OK" indicates an error. Typical error responses are "out of memory", "internal error", "error in input data" and so on. A non-OK response means there are no reliable solutions to look for.

Assuming the optimizer terminated without error, the next thing to ask is "why did it terminate?". The answer is provided by the termination code. "OK" means standard termination, reaching the numerical termination criteria at hand. Other termination codes, however, are not errors. They can indicate such reasons for stopping as "reached a given number of solutions", "reached the time limit", "could not make more progress (stall)", "user interruption from a callback function" and so on. In these cases there is usually a useful solution.

Having checked why the optimizer stopped the next question is "What does it think about my problem?". This is answered by the problem status, which, roughly, indicates if the problem was determined to be feasible or not. Problem status can also be unknown if no conclusion was reached.

Finally we can ask "How do I interpret the solution the optimizer returns?". This is indicated by the solution status. Besides a standard primal-dual feasible solution the solver can for instance return a certificate of primal or dual infeasibility. For integer problems solution status will indicate whether the solution is optimal or just feasible. Again, "unknown" is always a possible answer.

If we are happy with problem and solution status it is now time to fetch the results.

Some users will next want to analyze the solution on their own to check that it is numerically satisfying in that violations and gaps are sufficiently small. The last word on whether the solution is satisfactory always belongs to the user; the answer may be good enough even when the problem status is unknown and/or the termination code indicates "stall".

Finally some users may want to know in greater detail how the optimization process proceeded: data such as time taken by the optimizer, number of iterations, solution quality, etc. which normally appear in the log are available as information items.

This completes the discussion mostly from the perspective of the Optimizer API. In other interfaces implementation details will differ, for instance an exception will be thrown in place of a non-OK response code. Details and an example can be found in the "Accessing a solution" section of your API's tutorial and the corresponding response.* example.

Monday, February 18, 2019

Porting Fusion for MATLAB code

Fusion for MATLAB is discontinued starting with MOSEK version 9. If you wish to continue using MOSEK in MATLAB through the Fusion interface, this document describes how to adapt existing code to use the standard Fusion for Java.

First, use mosek.jar instead of mosekmatlab.jar in the Java path in MATLAB, for instance:
Next, every time you explicitly index into an object of a class from mosek.fusion use 0-based indexing, standard for Java, instead of 1-based indexing familiar from MATLAB. This applies to classes such as Variable, Expression, Matrix and to all operations such as index, pick, slice. For example, if you define a variable  then its individual entries are
whereas in the old interface they would be indexed from 1 through 4.

This applies only to indexing Java objects from mosek.fusion. Nothing changes with regard to MATLAB arrays, so operations which don't use explicit indexes are not affected, and whenever the input/output of a method is a MATLAB array, it will be 1-based indexed as always. So the following piece of code works both in the old and new regime:

Finally note that Java 1.8+ is required.

Wednesday, January 30, 2019

Planned end of support for version 7

In connection with the recent beta release of MOSEK 9, we will be stopping support for version 7 which was released in 2013. According to our policy of supporting version $n$ at least two years after the release of version $n+1$, we plan to:

  • end support for version 7 on 31-jun-2019 (ie. over three years after the release of version 8)
  • support version 8 (current stable release) until at least 15-jan-2021 (and potentially longer)

Friday, January 25, 2019

DTU CEE Energy Systems Summer School 2019

DTU CEE Summer School 2019 "Data-Driven Analytics and Optimization for Energy Systems" takes place June 16-21 2019 at DTU in Copenhagen. This is the 4th edition in a very successful series of workshops about optimization in power systems. MOSEK sponsors the school and two MOSEK scholarships will be awarded to outstanding student applicants.

Program, speakers, deadlines and more can be found at: