1% 2% Toolkit of useful utilities for CLP(BNR) 3% 4/* The MIT License (MIT) 5 * 6 * Copyright (c) 2022-2024 Rick Workman 7 * 8 * Permission is hereby granted, free of charge, to any person obtaining a copy 9 * of this software and associated documentation files (the "Software"), to deal 10 * in the Software without restriction, including without limitation the rights 11 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 12 * copies of the Software, and to permit persons to whom the Software is 13 * furnished to do so, subject to the following conditions: 14 * 15 * The above copyright notice and this permission notice shall be included in all 16 * copies or substantial portions of the Software. 17 * 18 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 19 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 20 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 21 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 22 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 23 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 24 * SOFTWARE. 25 */ 26:- module(clpBNR_toolkit, % SWI module declaration 27 [ 28 iterate_until/3, % general purpose iterator 29 mid_split_one/1, % contractor to split largest interval at midpoint 30 mid_split/1, % contractor to split an interval at midpoint 31 taylor_contractor/2, % build cf_contractor based on Taylor expansion 32 taylor_merged_contractor/2, % build merged Taylor cf_contractor from list of equations 33 cf_contractor/2, % execute cf_contractor 34 cf_solve/1, cf_solve/2, % a solve predicate for centre form contractors 35 36 lin_minimum/3, % find minimum of linear problem using library(simplex) 37 lin_maximum/3, % find maximum of linear problem using library(simplex) 38 lin_minimize/3, % lin_minimum/3 plus bind vars to solution minimizers 39 lin_maximize/3, % lin_maximum/3 plus bind vars to solution maximizers 40 41 local_minima/1, % apply KT constraints for objective function expression (OFE) 42 local_maxima/1, % semantically equivalent to local_minima/1 43 local_minima/2, % apply KT constraints for minima with constraints 44 local_maxima/2 % apply KT constraints for maxima with constraints 45 ]).
56:- use_module(library(apply),[maplist/3]). 57:- use_module(library(clpBNR)). 58:- use_module(library(simplex)). 59 60% messages for noisy failure 61fail_msg_(FString,Args) :- 62 debug(clpBNR,FString,Args), fail. 63 64:- set_prolog_flag(optimise,true). % for arithmetic, this module only
small/2 and Goal midsplit/1 :
?- X::real(-1,1),iterate_until(10,small(X,0),mid_split(X)),format("X = ~w\n",X),fail;true.
X = _6288{real(-1,-1r2)}
X = _6288{real(-1r2,0)}
X = _6288{real(0,1r2)}
X = _6288{real(1r2,1)}
true.
The specific intended use case is to provide an iterator for meta-contractors such as the centre-form contractor such as midsplit/1 (example above) or as constructed by taylor_contractor/2 as in:
?- X::real,taylor_contractor({X**4-4*X**3+4*X**2-4*X+3==0},T),
iterate_until(50,small(X),(T,mid_split_one([X]))),format("X = ~w\n",X),fail;true.
X = _150{real(0.999999999926943,1.00000000007306)}
X = _150{real(2.999999999484828,3.0000000005152105)}
true.
(Aside: For some problems, solving with Taylor contractors can be a faster and more precise alternative to clpBNR:solve/1.)
*/
88% 89% General purpose iterator: execute Goal a maximum of N times or until Test succeeds 90% 91iterate_until(N,Test,Goal) :- N>0, !, 92 , 93 N1 is N-1, 94 ( 95 -> true 96 ; iterate_until(N1,Test,Goal) 97 ). 98iterate_until(_N,_,_). % non-positive N --> exit 99 100sandbox:safe_meta(clpBNR_toolkit:iterate_until(_N,Test,Goal), [Test, Goal]).
mid_split for details of interval splitting for this predicate.
109mid_split_one(Xs) :- 110 select_split(Xs,X), % select largest interval with largest width 111 mid_split(X). % split it
mid_split(X) :-
M is midpoint(X),
({X=<M} ; {M=<X}).
Note that mid_split succeeds if X is a number, but doesn't do anything.
Use clpBNR:small as a pre-test to avoid splitting intervals which are already small enough.
127mid_split(X) :- 128 (number(X) % optimise number case 129 -> true 130 ; (small(X) 131 -> true 132 ; midpoint(X,M), % fails if not an interval 133 ({X=<M} ; {M=<X}) % possible choicepoint 134 ) 135 ). 136% 137% select interval with largest width 138% 139select_split([X],X) :- !. % select last remaining element 140select_split([X1,X2|Xs],X) :- % compare widths and discard one interval 141 delta(X1,D1), 142 delta(X2,D2), 143 (D1 >= D2 144 -> select_split([X1|Xs],X) 145 ; select_split([X2|Xs],X) 146 ).
taylor_contractor. In normal usage, a direct call to cf_contractor does appear; instead use cf_contractor or in a Goal for iterate_until/3.
155% 156% centred form contractor 157% 158% Bind the values of As to the midpoints of Xs. To support repetitive application 159% of the contractor (required by the iterator), the contractor should not permanently 160% bind anything so findall/3 will be used to achieve this "forward checking" 161% (as suggested in [CLIP]). After the call to findall, the bounds of the resulting list 162% of narrowed domains (XDs) are then applied to Xs. 163% 164% This contractor can be used with any "centred form", e.g., Newton or Krawczyk, since it 165% only depends on intervals and their midpoints, hence its name `cf_contractor`. The 166% details which distinguish the variety of centred form are built into the variables' 167% constraints. 168% 169cf_contractor(Xs,As) :- 170 findall(Ds,(maplist(bind_to_midpoint,Xs,As),maplist(cf_domain,Xs,Ds)),[XDs]), 171 maplist(set_domain,Xs,XDs). 172 173bind_to_midpoint(X,A) :- A is float(midpoint(X)). 174 175cf_domain(X,D) :- 176 number(X) -> D = X ; domain(X,D). % in case X narrowed to a point 177 178set_domain(X,D) :- 179 number(D) -> X = D ; X::D.
clpBNR_default_precision); otherwise fails.
This is done by using iterate_until/3 limited to a count determined by the flag clpBNR_iteration_limit. Examples:
?- X::real, taylor_contractor({X**4-4*X**3+4*X**2-4*X+3==0},T), cf_solve(T).
T = cf_contractor([X], [_A]),
X:: 1.000000000...,
_A::real(-1.0Inf, 1.0Inf) ;
T = cf_contractor([X], [_A]),
X:: 3.00000000...,
_A::real(-1.0Inf, 1.0Inf) ;
false.
?- taylor_contractor({2*X1+5*X1**3+1==X2*(1+X2), 2*X2+5*X2**3+1==X1*(1+X1)},T), cf_solve(T).
T = cf_contractor([X2, X1], [_A, _B]),
X1:: -0.42730462...,
X2:: -0.42730462...,
_B::real(-1.0Inf, 1.0Inf),
_A::real(-1.0Inf, 1.0Inf) ;
false.
211cf_solve(T) :- 212 current_prolog_flag(clpBNR_default_precision,P), 213 cf_solve(T,P). 214cf_solve(cf_contractor(Xs,As),P) :- 215 current_prolog_flag(clpBNR_iteration_limit,L), 216 Count is L div 10, % heuristic - primitive iteration limit/10 217 cf_iterate_(Count,Xs,As,P). 218 219cf_iterate_(Count,Xs,As,P) :- 220 Count > 0, 221 \+ small(Xs,P), % at least one var not narrow enough 222 !, 223 cf_contractor(Xs,As), % execute contractor 224 select_split(Xs,X), % select widest 225 (small(X,P) % still wide enough to split? 226 -> true % no, we're done 227 ; mid_split(X), % yes, split it 228 Count1 is Count-1, 229 cf_iterate_(Count1,Xs,As,P) % and iterate 230 ). 231cf_iterate_(_,_,_,_). % done (Count=<0 or all small Xs)
== or =:=) constraints Constraints; otherwise fails. Example:
?- taylor_contractor({X**4-4*X**3+4*X**2-4*X+3==0},T).
T = cf_contractor([X], [_A]),
X::real(-1.509169756145379, 4.18727500493995),
_A::real(-1.0Inf, 1.0Inf).
Use the contractor with cf_solve to search for solutions, as in:
?- X::real,taylor_contractor({X**4-4*X**3+4*X**2-4*X+3==0},T), cf_solve(T).
T = cf_contractor([X], [_A]),
X:: 1.000000000...,
_A::real(-1.0Inf, 1.0Inf) ;
T = cf_contractor([X], [_A]),
X:: 3.00000000...,
_A::real(-1.0Inf, 1.0Inf) ;
false.
Multiple equality constraints are supported, as in this example of the Broyden banded problem (N=2):
?- taylor_contractor({2*X1+5*X1**3+1==X2*(1+X2), 2*X2+5*X2**3+1==X1*(1+X1)},T), cf_solve(T).
T = cf_contractor([X2, X1], [_A, _B]),
X1:: -0.42730462...,
X2:: -0.42730462...,
_B::real(-1.0Inf, 1.0Inf),
_A::real(-1.0Inf, 1.0Inf) ;
false.
Centre form contractors can converge faster than the general purpose builtin fixed point iteration provided by solve/1.
268% 269% build a cf_contractor for a multivariate expression based on Taylor expansion 270% 271taylor_contractor({E1=:=E2},CF) :- 272 taylor_contractor({E1==E2},CF). 273taylor_contractor({E1==E2},cf_contractor(Xs,As)) :- 274 Exp=E1-E2, 275 term_variables(Exp,Xs), % original arguments, bound to TXs on call 276 make_EQ_(Exp,TEQ), % original constraint with arguments 277 % build constraint list starting with Z's and ending with TEQ and DEQ () 278 T::real(0,1), 279 make_As_and_Zs_(Xs,T,As,Zs,Cs,[TEQ,DEQ]), % T per Z 280 % now build Taylor constraint, DEQ 281 copy_term_nat(Exp,AExp), % copy of original constraint with As 282 term_variables(AExp,As), 283 sum_diffs(Xs, As, Zs, Zs, Exp, AExp, DEQ), % add on D(Z)'s' 284 % make any vars in original equation and contractor arguments finite real intervals 285 !, 286 Xs::real, % all vars are intervals 287 {Cs}. % apply constraints 288taylor_contractor({Es},CF) :- 289 taylor_merged_contractor({Es},CF), % list or sequence 290 !. 291taylor_contractor(Eq,_) :- 292 fail_msg_('Invalid constraint for Taylor contractor: ~w',[Eq]). 293 294make_As_and_Zs_([],_,[],[],Tail,Tail). 295make_As_and_Zs_([X|Xs],T,[A|As],[Z|Zs],[Z==A+T*(X-A)|CZs],Tail) :- 296 make_As_and_Zs_(Xs,T,As,Zs,CZs,Tail). 297 298sum_diffs([], [], [], _AllZs, _Exp, ExpIn, EQ) :- make_EQ_(ExpIn,EQ). 299sum_diffs([X|Xs], [A|As], [Z|Zs], AllZs, Exp, AExp, DEQ) :- 300 copy_term_nat(Exp,NExp), % copy expression and replace Xs by Zs 301 term_variables(NExp,AllZs), 302 partial_derivative(NExp,Z,DZ), % differentiate wrt. Z and add to generated expression 303 sum_diffs(Xs, As, Zs, AllZs, Exp, AExp+DZ*(X-A), DEQ). 304 305% map expression Exp to an equation equivalent to Exp==0 with numeric RHS 306make_EQ_(Exp,LHS==RHS) :- % turn expression into equation equivalent to Exp==0. 307 make_EQ_(Exp,LHS,RHS). 308 309make_EQ_(E,E,0) :- var(E), !. 310make_EQ_(X+Y,X,SY) :- number(Y), !, SY is -Y. 311make_EQ_(X-Y,X,Y) :- number(Y), !. 312make_EQ_(X+Y,Y,SX) :- number(X), !, SX is -X. 313make_EQ_(X-Y,SY,SX) :- number(X), !, SX is -X, negate_sum_(Y,SY). 314make_EQ_(X+Y,LHS+Y,RHS) :- !, make_EQ_(X,LHS,RHS). 315make_EQ_(X-Y,LHS-Y,RHS) :- !, make_EQ_(X,LHS,RHS). 316make_EQ_(E,E,0). % default (non +/- subexpression) 317 318negate_sum_(Y,-Y) :- var(Y), !. 319negate_sum_(X+Y,NX-Y) :- !, negate_sum_(X,NX). 320negate_sum_(X-Y,NX+Y) :- !, negate_sum_(X,NX). 321negate_sum_(E,-E).
== or =:=) constraint in Constraints; otherwise fails.
330% 331% build a cf_contractor by merging a list of cf_contractor's 332% 333taylor_merged_contractor({Es},T) :- 334 cf_list(Es,Ts), 335 cf_merge(Ts,T). 336 337cf_list([],[]) :- !. 338cf_list([C|Cs],[CF|CFs]) :- !, 339 cf_list(C, CF), 340 cf_list(Cs,CFs). 341cf_list((C,Cs),[CF|CFs]) :- !, 342 cf_list(C, CF), 343 cf_list(Cs,CFs). 344cf_list(C,CF) :- 345 taylor_contractor({C},CF). 346 347cf_merge(CFs,CF) :- cf_merge(CFs,cf_contractor([],[]),CF). 348 349cf_merge([],CF,CF). 350cf_merge([CF|CFs],CFIn,CFOut) :- 351 cf_merge(CF,CFIn,CFNxt), 352 cf_merge(CFs,CFNxt,CFOut). 353cf_merge(cf_contractor(Xs,As),cf_contractor(XsIn,AsIn),cf_contractor(XsOut,AsOut)) :- 354 cf_add(Xs,As,XsIn,AsIn,XsOut,AsOut). 355 356cf_add([],[],Xs,As,Xs,As). 357cf_add([X|Xs],[A|As],XsIn,AsIn,XsOut,AsOut) :- 358 var_existing(XsIn,AsIn,X,A), !, 359 cf_add(Xs,As,XsIn,AsIn,XsOut,AsOut). 360cf_add([X|Xs],[A|As],XsIn,AsIn,XsOut,AsOut) :- 361 cf_add(Xs,As,[X|XsIn],[A|AsIn],XsOut,AsOut). 362 363var_existing([Xex|Xs],[Aex|As], X,A) :- Xex==X -> Aex=A ; var_existing(Xs,As,X,A).
X*C (or C*X) are permitted since the actual computation is done using library(simplex). Narrowing of minimizers (variables in ObjF) is limited to that constrained by the Min result to accomodate multiple sets of minimizers. (See lin_minimize/3 to use minimizers used to derive Min.) A solution generator, e.g., clpBNR:solve/1 can be used to search for alternative sets of minimizers. "Universal Mines" example from the User Guide:
?- [M_Idays,M_IIdays,M_IIIdays]::integer(0,7),
lin_minimum(20*M_Idays+22*M_IIdays+18*M_IIIdays,
{4*M_Idays+6*M_IIdays+M_IIIdays>=54,4*M_Idays+4*M_IIdays+6*M_IIIdays>=65}, Min).
Min = 284,
M_Idays::integer(2, 7),
M_IIdays::integer(4, 7),
M_IIIdays::integer(2, 7).
?- [M_Idays,M_IIdays,M_IIIdays]::integer(0,7),
lin_minimum(20*M_Idays+22*M_IIdays+18*M_IIIdays,
{4*M_Idays+6*M_IIdays+M_IIIdays>=54,4*M_Idays+4*M_IIdays+6*M_IIIdays>=65}, Min),
solve([M_Idays,M_IIdays,M_IIIdays]).
M_Idays = 2,
M_IIdays = 7,
M_IIIdays = 5,
Min = 284 ;
false.
For linear systems, lin_minimum/3, lin_maximum/3 can be significantly faster than using the more general purpose clpBNR:global_minimum/3, clpBNR:global_maximum/3
lin_minimum/3 for finding global maxima.
400lin_minimum(ObjF,{Constraints},MinValue) :- 401 lin_minimum_(ObjF,{Constraints},MinValue,false). 402 403lin_maximum(ObjF,{Constraints},MinValue) :- 404 lin_maximum_(ObjF,{Constraints},MinValue,false).
lin_minimum/3 except variables in ObjF will be narrowed to the values used in calculating the final value of Min. Any other sets of minimizers corresponding to Min are removed from the solution space. "Universal Mines" example from the User Guide:
?- [M_Idays,M_IIdays,M_IIIdays]::integer(0,7),
lin_minimize(20*M_Idays+22*M_IIdays+18*M_IIIdays,
{4*M_Idays+6*M_IIdays+M_IIIdays>=54,4*M_Idays+4*M_IIdays+6*M_IIIdays>=65}, Min).
M_Idays = 2,
M_IIdays = 7,
M_IIIdays = 5,
Min = 284.
lin_maximum/3 except variables in ObjF will be narrowed to the values used in calculating the final value of Max. Any other sets of minimizers corresponding to Min are removed from the solution space.
427lin_minimize(ObjF,{Constraints},MinValue) :- 428 lin_minimum_(ObjF,{Constraints},MinValue,true). 429 430lin_maximize(ObjF,{Constraints},MinValue) :- 431 lin_maximum_(ObjF,{Constraints},MinValue,true). 432 433 434lin_minimum_(ObjF,{Constraints},MinValue,BindVars) :- 435 init_simplex_(ObjF,Constraints,Objective,S0,Vs), 436 (minimize(Objective,S0,S) 437 -> objective(S,MinValue), {ObjF == MinValue}, 438 (BindVars == true 439 -> bind_vars_(Vs,S) 440 ; remove_names_(Vs), 441 {Constraints} % apply constraints 442 ) 443 ; fail_msg_('Failed to minimize: ~w',[ObjF]) 444 ). 445 446lin_maximum_(ObjF,{Constraints},MaxValue,BindVars) :- 447 init_simplex_(ObjF,Constraints,Objective,S0,Vs), 448 (maximize(Objective,S0,S) 449 -> objective(S,MaxValue), {ObjF == MaxValue}, 450 (BindVars == true 451 -> bind_vars_(Vs,S) 452 ; remove_names_(Vs), 453 {Constraints} % apply constraints 454 ) 455 ; fail_msg_('Failed to maximize: ~w',[ObjF]) 456 ). 457 458init_simplex_(ObjF,Constraints,Objective,S,Vs) :- 459 gen_state(S0), 460 term_variables((ObjF,Constraints),Vs), 461 (Vs = [] 462 -> fail_msg_('No variables in expression to optimize: ~w',[ObjF]) 463 ; sim_constraints_(Constraints,S0,S1), 464 _::real(_,Max), % max value to constrain for simplex 465 constrain_ints_(Vs,Max,S1,S), 466 (map_simplex_(ObjF,T/T,Objective/[]) 467 -> true 468 ; fail_msg_('Illegal linear objective: ~w',[ObjF]) 469 ) 470 ). 471 472% use an attribute to associate a var with a unique simplex variable name 473simplex_var_(V,SV) :- var(V), 474 (get_attr(V,clpBNR_toolkit,SV) 475 -> true 476 ; statistics(inferences,VID), SV = var(VID), put_attr(V,clpBNR_toolkit,SV) 477 ). 478 479% Name attribute removed before exit. 480remove_names_([]). 481remove_names_([V|Vs]) :- 482 del_attr(V,clpBNR_toolkit), 483 remove_names_(Vs). 484 485attr_unify_hook(var(_),_). % unification always does nothing and succeeds 486 487constrain_ints_([],_,S,S). 488constrain_ints_([V|Vs],Max,Sin,Sout) :- 489 % Note: library(simplex) currently constrains all variables to be non-negative 490 simplex_var_(V,SV), % corresponding simplex variable name 491 % if not already an interval, make it one with finite non-negative value 492 (domain(V,D) -> true ; V::real(0,_), domain(V,D)), 493 (D == boolean -> Dom = integer(0,1); Dom = D), 494 Dom =.. [Type,L,H], 495 (Type == integer -> constraint(integral(SV),Sin,S1) ; S1 = Sin), 496 (L < 0 497 -> % apply non-negativity condition 498 ({V >= 0} -> L1 = 0 ; fail_msg_('Negative vars not supported by \'simplex\': ~w',[V])) 499 ; L1 = L 500 ), 501 % explicitly constrain any vars not (0,Max-eps) 502 (L1 > 0 -> constraint([SV] >= L1,S1,S2) ; S2 = S1), % L1 not negative 503 (H < Max -> constraint([SV] =< H,S2,SNxt) ; SNxt = S2), 504 constrain_ints_(Vs,Max,SNxt,Sout). 505 506bind_vars_([],_). 507bind_vars_([V|Vs],S) :- 508 % Note: skip anything nonvar (already bound due to active constraints) 509 (simplex_var_(V,SV) -> variable_value(S,SV,V) ; true), 510 bind_vars_(Vs,S). 511 512% clpBNR constraints have already been applied so worst errors have been detected 513sim_constraints_([],S,S) :- !. 514sim_constraints_([C|Cs],Sin,Sout) :- !, 515 sim_constraints_(C, Sin,Snxt), 516 sim_constraints_(Cs,Snxt,Sout). 517sim_constraints_((C,Cs),Sin,Sout) :- !, 518 sim_constraints_(C, Sin,Snxt), 519 sim_constraints_(Cs,Snxt,Sout). 520sim_constraints_(C,Sin,Sout) :- 521 sim_constraint_(C,SC), 522 constraint(SC,Sin,Sout). % from simplex 523 524sim_constraint_(C,SC) :- 525 C=..[Op,LHS,RHS], % decompose 526 constraint_op(Op,COp), % acceptable to simplex 527 number(RHS), RHS >= 0, % requirement of simplex 528 map_simplex_(LHS,T/T,Sim/[]), % map to simplex list of terms 529 !, 530 SC=..[COp,Sim,RHS]. % recompose 531sim_constraint_(C,_) :- 532 fail_msg_('Illegal linear constraint: ~w',[C]). 533 534map_simplex_(T,CT/[S|Tail],CT/Tail) :- 535 map_simplex_term_(T,S), 536 !. 537map_simplex_(A+B,Cin,Cout) :- !, 538 map_simplex_(A,Cin,Cnxt), 539 map_simplex_(B,Cnxt,Cout). 540map_simplex_(A-B,Cin,Cout) :- !, 541 map_simplex_(A,Cin,Cnxt), 542 map_simplex_(-B,Cnxt,Cout). 543 544map_simplex_term_(V,1*SV) :- simplex_var_(V,SV), !. 545map_simplex_term_(-T,NN*V) :- !, 546 map_simplex_term_(T,N*V), 547 NN is -N. 548map_simplex_term_(N*V,N*SV) :- number(N), simplex_var_(V,SV), !. 549map_simplex_term_(V*N,N*SV) :- number(N), simplex_var_(V,SV). 550 551constraint_op(==,=). 552constraint_op(=<,=<). 553constraint_op(>=,>=).
local_minima should be executed prior to a call to clpBNR:global_minimum using the same objective function, e.g.,
?- X::real(0,10), OF=X**3-6*X**2+9*X+6, local_minima(OF), global_minimum(OF,Z). OF = X**3-6*X**2+9*X+6, X:: 3.00000000000000..., Z:: 6.000000000000... .
Using any local optima predicate can significantly improve performance compared to searching for global optima (clpBNR:global_*) without local constraints.
local_maxima should be executed prior to a call to clpBNR:global_maximum using the same objective function, e.g.,
?- X::real(0,10), OF=X**3-6*X**2+9*X+6, local_maxima(OF), global_maximum(OF,Z). OF = X**3-6*X**2+9*X+6, X:: 1.000000000000000..., Z:: 10.0000000000000... .
581% 582% local_minima/1, % apply KT constraints for objective function expression (OFE) 583% local_maxima/1, % semantically equivalent to local_minima/1 584% 585local_minima(ObjExp) :- 586 local_optima_(min,ObjExp,[]). 587 588local_maxima(ObjExp) :- 589 local_optima_(max,ObjExp,[]).
local_minima should be executed prior to a call to clpBNR:global_minimum using the same objective function, e.g.,
?- [X1,X2]::real, OF=X1**4*exp(-0.01*(X1*X2)**2),
local_minima(OF,{2*X1**2+X2**2==10}), global_minimum(OF,Z), solve([X1,X2]).
OF = X1**4*exp(-0.01*(X1*X2)**2),
X1::real(-1.703183936003284e-108, 1.703183936003284e-108),
X2:: -3.16227766016838...,
Z:: 0.0000000000000000... ;
OF = X1**4*exp(-0.01*(X1*X2)**2),
X1::real(-1.703183936003284e-108, 1.703183936003284e-108),
X2:: 3.16227766016838...,
Z:: 0.0000000000000000... .
local_maxima should be executed prior to a call to clpBNR:global_maximum using the same objective function, e.g.,
?- [X1,X2]::real,OF=X1**4*exp(-0.01*(X1*X2)**2),
local_maxima(OF,{2*X1**2+X2**2==10}), global_maximum(OF,Z),solve([X1,X2]).
OF = X1**4*exp(-0.01*(X1*X2)**2),
X1:: -2.23606797749979...,
X2:: 0.0000000000000000...,
Z:: 25.0000000000000... ;
OF = X1**4*exp(-0.01*(X1*X2)**2),
X1:: 2.23606797749979...,
X2:: 0.0000000000000000...,
Z:: 25.0000000000000... .
628% 629% local_minima/2, % apply KT constraints for minima with constraints 630% local_maxima/2 % apply KT constraints for maxima with constraints 631% 632local_minima(ObjExp,{Constraints}) :- 633 local_optima_(min,ObjExp,Constraints). 634 635local_maxima(ObjExp,{Constraints}) :- 636 local_optima_(max,ObjExp,Constraints). 637 638 639local_optima_(MinMax,ObjF,Constraints) :- 640 local_optima_(MinMax,ObjF,Constraints,Cs), % generate constraints 641 {Cs}. % then apply 642 643local_optima_(MinMax,ObjF,Constraints,[Constraints,GCs,LCs]) :- 644 lagrangian_(Constraints,MinMax,ObjF,LObjF,LCs), 645 term_variables((Constraints,ObjF),Vs), 646 gradient_constraints_(Vs,GCs,LObjF). 647 648gradient_constraints_([],[],_Exp). 649gradient_constraints_([X|Xs],[C|Cs],Exp) :- 650 partial_derivative(Exp,X,D), 651 (number(D) -> C=[] ; C=(D==0)), % no constraint if PD is a constant 652 gradient_constraints_(Xs,Cs,Exp). 653 654% for each constraint add to Lagrangian expression with auxiliary KKT constraints 655lagrangian_(C,MinMax,Exp,LExp, LC) :- nonvar(C), 656 kt_constraint_(C,CExp, LC), % generate langrange term with multiplier 657 lexp(MinMax,Exp,CExp,LExp), 658 !. 659lagrangian_([],_,L,L,[]). 660lagrangian_([C|Cs],MinMax,Exp,LExp,[LC|LCs]) :- 661 lagrangian_(C, MinMax,Exp,NExp,LC), 662 lagrangian_(Cs,MinMax,NExp,LExp,LCs). 663lagrangian_((C,Cs),MinMax,Exp,LExp,[LC|LCs]) :- 664 lagrangian_(C,MinMax,Exp,NExp,LC), 665 lagrangian_(Cs,MinMax,NExp,LExp,LCs). 666 667lexp(min,Exp,CExp,Exp+CExp). 668lexp(max,Exp,CExp,Exp-CExp). 669 670kt_constraint_(LHS == RHS, M*(LHS-RHS), []) :- 671 M::real. % finite multiplier only 672kt_constraint_(LHS =< RHS, MGx, MGx==0) :- 673 MGx = M*(LHS-RHS), M::real(0,_). % positive multiplier and KKT non-negativity condition 674kt_constraint_(LHS >= RHS, Exp, LC) :- 675 kt_constraint_(RHS =< LHS, Exp, LC). % >= convert to =<
clpBNR_toolkit: Toolkit of various utilities used for solving problems with clpBNR
CLP(BNR) (
library(clpBNR))is a CLP over the domain of real numbers extended with ±∞. This module contains a number of useful utilities for specific problem domains like the optimization of linear systems, enforcing local optima conditions, and constructing centre form contractors to improve performance (e.g., Taylor extensions of constraints). For more detailed discussion, see A Guide to CLP(BNR) (HTML version included with this pack in directorydocs/).Documentation for exported predicates follows. The "custom" types include:
clpBNRattribute