Documentation for global_maximum/2 etc?

[JeanChristopheRohner]

Thanks for porting clp(BNR)! I read that the documentation is work in progress at this point, but is there some place where i can find information about how to use global_maximum/2 and related predicates?

[ridgeworks]

Thanks for the query. As it happens I just finished a draft on global optimization for the Guide, but I'm not ready to "publish" it yet.

I've attached a section of the document source which should help you get started. It's in a markup language similar to Markdown so I don't think you should have any problems understanding it. (You may have to soft wrap the text in your editor.)

Hope this keeps you going until the Guide gets updated, but let me know if you have any specific questions.

On Jul 28, 2020, at 5:00 AM, JeanChristopheRohner notifications@github.com wrote:

Thanks for porting clp(BNR)! I read that the documentation is work in progress at this point, but is there some place where i can find information about how to use global_maximum/2 and related predicates?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/ridgeworks/clpBNR_pl/issues/2, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABXFJQPSGYKRZFVI6VNFPV3R52HS7ANCNFSM4PKJLM6A.

Global Optimization

An important class of problems in engineering design and operations research is the optimal allocation of "resources", usually referred to as global optimization. Specifically it aims to find the maximum or minimum value of a function of one or more variables and to find the values of the variables where this optima is reached. Techniques for solving local optimization problems, i.e., solutions in the neighbourhood of a given point, have been well-researched but techniques for global optimization seem to be less well developed. One of the reasons for this is the difficulty of representing global information, but intervals can be used to address this issue. The evaluation of a function of intervals using interval arithmetic is an interval containing all the values of the function given the interval values of the arguments. Then using search techniques similar to solving for the roots of equations, a minimum/maximum value of the function, along with it's minimizers (the values of the arguments producing the minimum/maximum) can be found.

Such searches can be computationally expensive, so it is often advantageous to combine it with local optimization techniques. These take the form of additional (redundant) constraints which can accelerate the searching process. In addition many global optimization problems include additional problem constraints which put further conditions on the set of solutions, so the combination of Constraint Logic Programming and Relational Interval Arithmetic seems to be a good fit for tackling this class of problems.

Material in this section draws heavily on two sources: for theory, see [New Computer Methods for Global Optimization] (Ratschek and Rokne, 2007). The examples are largely from [Schaum's Outlines, Operations Research] (Bronson and Naadimuthu,1997).

Unconstrained Optimization

CLP(BNR) provides a pair of general purpose predicates for tackling global optimization problems: global_minimum and global_maximum. The simplest form of these predicates finds the global min/max (the second argument, Z, in the examples below) for the objective function specified by the first argument. From "Shaum's", example 10.1: eg ?- X::real(0,20), global_maximum(X(5pi-X),Z). X:: 7.8..., Z:: 61.6850... . The predicates provided use the Moore-Skelboe algorithm (see Ratschek&Rokne) to find global optima to the precision specified by the clpBNR_default_precision flag (default 6 digits). They do not provide the minimisers - indeed there may be more than one minimiser value. In the example above some narrowing of X has occurred but, in general, one of the solve predicates can be used to separate and sharpen minimizer values. Since there really aren't any "point solutions" to the objective function (unlike like finding roots of an equation discussed earlier), splitsolve is usually better suited to this problem. eg ?- X::real(0,20), global_maximum(X(5pi-X),Z), splitsolve(X). X:: 7.84954..., Z:: 61.6850... . And, as we saw earlier, there may be answers for the same solution due to precision limits on the splitting process: eg ?- X::real(0,20), global_maximum(X(5pi-X),Z),splitsolve(X). X:: 7.84954..., Z:: 61.6850... ; X:: 7.84962..., Z:: 61.6850... ; X:: 7.8496..., Z:: 61.6850... ; X:: 7.84966..., Z:: 61.6850... .

global_minimum and global_maximum take an optional third argument: a precision value which overrides the current value of the clpBNR_default_precision flag. Since the underlying search is exhastive and can take some time dependent on the objective function and the domain(s) of the input, the precision argument can control the tradeoff between precision and execution time (as it does for other search predicates like solve and splitsolve): eg ?- X::real(0,20),global_maximum(X(5pi-X),Z,3). X::real(7.715468219981181, 7.992506304532757), Z:: 61.7... .

These predicates are general in nature; they can find optima on boundaries (where no local optima exist) and on functions of integer variables. Some examples of objective functions of more than one variable: eg ?- [X1,X2,X3]::real,global_maximum(- (2X1-5)2-(X2-3)2-(5X3-2)**2,Z). X1:: 2.500..., Z::real(-1.8040366808093222e-7, 0), X3:: 0.400..., X2:: 3.000... .

?- [X1,X2]::integer(-100,100),global_minimum((X1-sqrt(5))**2+(X2-pi)**2+10,Z).
X1 = 2,
X2 = 3,
Z:: 10.07577656955144... .

?- [X1,X2]::integer(-100,100),global_minimum((1-X1)**2+100*(X2-X1**2)**2,Z).
X1 = X2, X2 = 1,
Z = 0.

Constrained Optimization

Moving to CLP usually involves a paradigm shift from "generate and test" to to "apply test, then generate". This facilitates solving constrained optimization problems; just apply the constraints before calling the global optimization predicate. Some examples from "Shaum's" (chapter 12): eg ?- {X2+Y2==1},global_minimum(X+Y,Z). X:: -0.70..., Y:: -0.70..., Z:: -1.41421... .

?- {X**2+Y**2==1},global_minimum(X+Y,Z),splitsolve([X,Y]).
X:: -0.707107...,
Y:: -0.70710...,
Z:: -1.41421... .

?- [X1,X2,X3]::real(0,2*pi), {-X1*X2**3+X1**2*X3**2==5}, global_maximum(sin(X1*X2+X3),Z),splitsolve([X1,X2,X3]).
X1:: 1.49526...,
X2:: 0.049725...,
Z:: 0.99999951...,
X3:: 1.49546... .

?- [X1,X2,X3]::integer,{X1+2*X2+3*X3==4,X1*X3==19}, global_maximum(-X1**6*X2**2-X1**4*X3**2-1,Z),splitsolve([X1,X2,X3]).
X1 = 1,
X2 = -27,
X3 = 19,
Z = -1091.

?- [X1,X2]::real(_,_),{X1+X2=<2}, global_maximum(log(1+X1)+2*log(1+X2),Z),splitsolve([X1,X2]).
X1:: 0.333333...,
X2:: 1.666666...,
Z:: 2.249340... .

?- [X1,X2]::real(0,_),{X1+2*X2=<3,8*X1+5*X2>=10}, global_minimum((X1-2)**2+(X2-2)**2,Z),splitsolve([X1,X2]).
X1:: 1.3992345...,
X2:: 0.8003827...,
Z:: 1.800000... .

Most of these examples generate answers fairly quickly, but it should be noted that the global optima predicates provided by CLP(BNR) are exhaustive searches. It's not hard to come up with examples that take a long time and millions of narrowing operations. Constraining the search space and using the tunable precision flag can help to produce acceptable answers in a reasonable time frame. And mixing local optimization techniques with global searches can result in sharper answers in less time as we'll see later.

[JeanChristopheRohner]

Thanks a lot for your answer. To be able to optimise non-linear expressions is very useful (clpfd could do so with integers but not with reals, clpqr with reals but only on linear expressions).

[ridgeworks]

Global optimization section added to User Guide with v0.9.1 (12/08/20).

[erikkaplun]

I have noticed than using global_minimum and global_maximum within the same scope on the same expression, the expression needs to be stored twice, i.e. in separate variables, because the variable holding the expression gets overwritten with the computed minimum or maximum value. This results in slightly confusing code. Is this planned? Here's an example that uses a completely brainless expression to give the library a try:

expression(X, Y, Z) :-
    { Z == exp(1/2)*X - (cos(Z)/Z)**(1/3) - 1+log((Y+3/Z)/Z) }.

optimize_one(Max_X, Max_Y, Min_Z-Max_Z) :-
    X::real(0.1, Max_X),
    Y::real(0.1, Max_Y),
    %% global_minimum/global_maximum "eats up" its first argument,
    %% so we need two copies of the same expression
    expression(X, Y, Z1),
    expression(X, Y, Z2),
    clpBNR:global_minimum(Z1, Min_Z),
    clpBNR:global_maximum(Z2, Max_Z).

optimize_many(MinMax_Zs) :-
    Max_Xs = [0.5, 0.8, 1.1, 10, 150, 1222],
    Max_Ys = [0.3, 0.7, 0.9, 25, 250, 2500].
    concurrent_maplist(optimize_one, Max_Xs, Max_Ys, MinMax_Zs).

[ridgeworks]

I think that's the only logical behaviour. There isn't a single X,Y pair that results in both a global maximum and minimum, i.e., the minimizing values of X and Y cannot be unified with the maximizing values of X and Y. You should be able to use ";' to produce one or the other (on backtracking).

Hope I'm not missing something.

[erikkaplun]

OK, for the sake of clarity, I will eliminate the use of global_maximum and will only use global_minimum. Note that Z and Z2 on lines 3 and 4 are exactly equivalent.

minmax_a(Min_Z, Max_Z) :-
    [X, Y]::real(0.1, 0.9),
    { Z  == exp(1/2)*X - (cos(Z)/Z)**(1/3)   - 1+log((Y+3/Z)/Z) },
    { Z2 == exp(1/2)*X - (cos(Z2)/Z2)**(1/3) - 1+log((Y+3/Z2)/Z2) },
    clpBNR:global_minimum(Z,   Min_Z),
    clpBNR:global_minimum(-Z2, Max_Z).

minmax_b(Min_Z, Max_Z) :-
    [X, Y]::real(0.1, 0.9),
    { Z == exp(1/2)*X - (cos(Z)/Z)**(1/3) - 1+log((Y+3/Z)/Z) },
    clpBNR:global_minimum(Z,  Min_Z),
    clpBNR:global_minimum(-Z, Max_Z).

?- minmax_a(Min_Z, Max_Z).
Min_Z:: -29.9999...,
Max_Z:: -0.91568... .

?- minmax_b(Min_Z, Max_Z).
Min_Z:: -29.9999...,
Max_Z:: 29.9999... .

— why are minmax_a and minmax_b producing different results if a pure reading of the program would imply that they should be equivalent?

[ridgeworks]

I don't see anything that would imply that the two predicates are equivalent. Z and Z2 are not equivalent; they're two different logic variables subject to different constraints. The constraints happen to share two other variables, X and Y so in some sense, they are co-dependent, i.e., there must be common values for X and Y that permit both Z and -Z2 to be minimized.

But your example does raise another important issue which I tried to cover in the documentation, but it's easy to lose it. global_minimum and global_maximum work by performing an iterative search using the variables in the function being minimized. In your case that's just Z (or -Z) so that's not much to work with. Without actually trying it, you would be better to minimize exp(1/2)*X - (cos(Z)/Z)**(1/3) - 1+log((Y+3/Z)/Z), an expression which contains all the variables.

Finally you have picked a rather complex expression, in fact one that's recursive (Z is a function of itself); is that intentional? Furthermore, by default, Z will contain 0, which means cos(Z)/Z and (Y+3/Z)/Z are not finitely bounded. I have no idea how this function behaves and the builtin constraint engine obviously doesn't either:

?- [X, Y]::real(0.1, 0.9),{ Z  == exp(1/2)*X - (cos(Z)/Z)**(1/3)   - 1+log((Y+3/Z)/Z) }.
X::real(0.09999999999999999, 0.9000000000000001),
Y::real(0.09999999999999999, 0.9000000000000001),
Z::real(-1.0Inf, 1.0Inf).

I would suggest you start with something simpler so the problem of expression complexity doesn't interfere with understanding how global min/max works. Based on my own experience, this is anything but simple.

[erikkaplun]

Finally you have picked a rather complex expression, in fact one that's recursive (Z is a function of itself); is that intentional?

I just quickly came up with something based on the documented example to demonstrate to a friend that the library won't choke on non-trivial math. There is no meaning whatsoever in that "formula".

[erikkaplun]

Okay, how about this:

minimize_a(Min_Z, Max_Z) :-
    [X, Y]::real(0.1, 0.9),
    { Z == exp(1/2)*X - 1+log((Y+3/X)/X) },
    clpBNR:global_minimum(Z, Min_Z),
    [X2, Y2]::real(0.1, 0.9),
    { Z2 == exp(1/2)*X2 - 1+log((Y2+3/X2)/X2) },
    clpBNR:global_minimum(-Z2, Max_Z).

minimize_b(Min_Z, Max_Z) :-
    [X, Y]::real(0.1, 0.9),
    { Z == exp(1/2)*X - 1+log((Y+3/X)/X) },
    clpBNR:global_minimum(Z, Min_Z),
    clpBNR:global_minimum(-Z, Max_Z).

?- minimize_a(Min_Z, Max_Z).
Min_Z:: 1.82274...,
Max_Z:: -4.89821... .

?- minimize_b(Min_Z, Max_Z).
Min_Z:: 1.82274...,
Max_Z:: -1.82274... .

My question is still: why does calling global_minimum(VarContainingExpr, ComputedMinOfThatExpr) have a side effect on VarContainingExpr so that I cannot compute, say, global_maximum(VarContaininExpr, ComputedMaxOfThatExpr) later on without having to use a fresh version of the expression held by VarContainingExpr.

To be clear: the question is not about how minimize_a works, but about why minimize_b doesn't produce results that are identical to minimize_a — should I assume global_minimum is not pure?

If I change the order of clauses in minimize_b, the results will be different:

minimize_b(Min_Z, Max_Z) :-
    [X, Y]::real(0.1, 0.9),
    { Z == exp(1/2)*X - 1+log((Y+3/X)/X) },
    clpBNR:global_minimum(-Z, Max_Z),
    clpBNR:global_minimum(Z, Min_Z).

?- minimize_b(Min_Z, Max_Z).
Min_Z:: 4.89821...,
Max_Z:: -4.89821... .

[erikkaplun]

Or even a simpler question:

minimize :-
    [X, Y]::real(0.1, 0.9),
    { Z == exp(1/2)*X - 1+log((Y+3/X)/X) },
    writeln(Z),
    clpBNR:global_minimum(Z, Min_Z),
    writeln(Z),
    writeln(Min_Z).

then why does Z, or the 1st argument to global_minimum, get overwritten with what I requested to be stored in the 2nd argument of the predicate:

?- minimize.
_26890{real(0.5037642492953187,6.217190420527862)}
_116{real(1.8227410945780738,1.8227421422963976)}
_116{real(1.8227410945780738,1.8227421422963976)}
true.

[ridgeworks]

I think there's something fundamental you're missing. This:

{ Z == exp(1/2)*X - 1+log((Y+3/X)/X) }

establishes constraints between X,Y and Z. Any changes (narrowing of their values) to any one of them can potentially change the others.

Now when I try to find the minimum value of Z that changes X and Y, what global_minimum does is to find and constrain Z to it's minimum value. So it's going to narrow (not overwrite) in the process. Also note that Z and Min_Z are in fact the same variable (they've been unified in the process).

The question I think you really want to ask is:

?- [X, Y]::real(0.1, 0.9), global_minimum(exp(1/2)*X - 1+log((Y+3/X)/X), Min).
X::real(0.385062117866536, 0.9000000000000001),
Y::real(0.09999999999999999, 0.9000000000000001),
Min:: 1.82274... .

i.e., what's the minimum value of the expression given the ranges of X andY. Note that while X has narrowed slightly, global_minimum doesn't, in general, solve for the minimizers, i.e., what are the actual values of X and Y that produce the minimum value. For that (in general), you need a second solve step:

?- [X, Y]::real(0.1, 0.9), global_minimum(exp(1/2)*X - 1+log((Y+3/X)/X), Min),solve([X,Y]).
X:: 0.90000...,
Y:: 0.10000...,
Min:: 1.82274... .

Perhaps it's now clearer why asking for the both minimum and maximum values isn't usually what you want. What you're really asking is there an (X,Y) pair that both minimizes and maximizes the expression, which is seldom the case.

Finally, when you say "get overwritten with what I requested to be stored in the 2nd argument of the predicate" suggests you're not quite thinking in logic programming terminology. Variables are not something that represent storage or get overwritten. They're variables in the mathematical sense; they're either unknown (unbound variable) or a known value which does not change (immutable). (clpBNR intervals are just variables but are constrained to be a range of numbers and that range can only get narrower as more constraints are applied. )

Hope this helps.

[erikkaplun]

Hope this helps.

It does indeed! Thanks a bunch!