1#  ___________________________________________________________________________
2#
3#  Pyomo: Python Optimization Modeling Objects
4#  Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC
5#  Under the terms of Contract DE-NA0003525 with National Technology and
6#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
7#  rights in this software.
8#  This software is distributed under the 3-clause BSD License.
9#  ___________________________________________________________________________
10
11from math import pow
12from numpy import inf
13from numpy.linalg import norm
14from pyomo.contrib.trustregion.filterMethod import (
15    FilterElement, Filter)
16from pyomo.contrib.trustregion.helper import (cloneXYZ, packXYZ)
17from pyomo.contrib.trustregion.Logger import Logger
18from pyomo.contrib.trustregion.PyomoInterface import (
19    PyomoInterface, ROMType)
20
21def TRF(m, eflist, config):
22    """The main function of the Trust Region Filter algorithm
23
24    m is a PyomoModel containing ExternalFunction() objects Model
25    requirements: m is a nonlinear program, with exactly one active
26    objective function.
27
28    eflist is a list of ExternalFunction objects that should be
29    treated with the trust region
30
31    config is the persistent set of variables defined
32    in the ConfigBlock class object
33
34    Return:
35    model is solved, variables are at optimal solution or
36    other exit condition.  model is left in reformulated form, with
37    some new variables introduced in a block named "tR" TODO: reverse
38    the transformation.
39    """
40
41    logger = Logger()
42    filteR = Filter()
43    problem = PyomoInterface(m, eflist, config)
44    x, y, z = problem.getInitialValue()
45
46    iteration = -1
47
48    romParam, yr = problem.buildROM(x, config.sample_radius)
49    #y = yr
50    rebuildROM = False
51    xk, yk, zk = cloneXYZ(x, y, z)
52    chik = 1e8
53    thetak = norm(yr - yk,1)
54
55    objk = problem.evaluateObj(x, y, z)
56
57    while True:
58        if iteration >= 0:
59            logger.printIteration(iteration)
60            #print(xk)
61        # increment iteration counter
62        iteration = iteration + 1
63        if iteration > config.max_it:
64            print("EXIT: Maxmium iterations\n")
65            break
66
67        ######  Why is this here ###########
68        if iteration == 1:
69            config.sample_region = False
70        ################################
71
72        # Keep Sample Region within Trust Region
73        if config.trust_radius < config.sample_radius:
74            config.sample_radius = max(
75                config.sample_radius_adjust*config.trust_radius,
76                config.delta_min)
77            rebuildROM = True
78
79        #Generate a RM r_k (x) that is kappa-fully linear on sigma k
80        if(rebuildROM):
81            #TODO: Ask Jonathan what variable 1e-3 should be
82            if config.trust_radius < 1e-3:
83                problem.romtype = ROMType.linear
84            else:
85                problem.romtype = config.reduced_model_type
86
87            romParam, yr = problem.buildROM(x, config.sample_radius)
88            #print(romParam)
89            #print(config.sample_radius)
90
91
92
93        # Criticality Check
94        if iteration > 0:
95            flag, chik = problem.criticalityCheck(x, y, z, romParam)
96            if (not flag):
97                raise Exception("Criticality Check fails!\n")
98
99        # Save the iteration information to the logger
100        logger.newIter(iteration,xk,yk,zk,thetak,objk,chik,
101                       config.print_variables)
102
103        # Check for Termination
104        if (thetak < config.ep_i and
105            chik < config.ep_chi and
106            config.sample_radius < config.ep_delta):
107            print("EXIT: OPTIMAL SOLUTION FOUND")
108            break
109
110        # If trust region very small and no progress is being made,
111        # terminate. The following condition must hold for two
112        # consecutive iterations.
113        if (config.trust_radius <= config.delta_min and thetak < config.ep_i):
114            if subopt_flag:
115                print("EXIT: FEASIBLE SOLUTION FOUND ")
116                break
117            else:
118                subopt_flag = True
119        else:
120            # This condition holds for iteration 0, which will declare
121            # the boolean subopt_flag
122            subopt_flag = False
123
124        # New criticality phase
125        if not config.sample_region:
126            config.sample_radius = config.trust_radius/2.0
127            if config.sample_radius > chik * config.criticality_check:
128                config.sample_radius = config.sample_radius/10.0
129            config.trust_radius = config.sample_radius*2
130        else:
131            config.sample_radius = max(min(config.sample_radius,
132                                   chik*config.criticality_check),
133                               config.delta_min)
134
135        logger.setCurIter(trustRadius=config.trust_radius,
136                          sampleRadius=config.sample_radius)
137
138        # Compatibility Check (Definition 2)
139        # radius=max(kappa_delta*config.trust_radius*min(1,kappa_mu*config.trust_radius**mu),
140        #            delta_min)
141        radius = max(config.kappa_delta*config.trust_radius *
142                     min(1,
143                         config.kappa_mu*pow(config.trust_radius,config.mu)),
144                     config.delta_min)
145
146        try:
147            flag, obj = problem.compatibilityCheck(
148                x, y, z, xk, yk, zk, romParam, radius,
149                config.compatibility_penalty)
150        except:
151            print("Compatibility check failed, unknown error")
152            raise
153
154        if not flag:
155            raise Exception("Compatibility check fails!\n")
156
157
158        theNorm = norm(x - xk, 2)**2 + norm(z - zk, 2)**2
159        if (obj - config.compatibility_penalty*theNorm >
160            config.ep_compatibility):
161            # Restoration stepNorm
162            yr = problem.evaluateDx(x)
163            theta = norm(yr - y, 1)
164
165            logger.iterlog.restoration = True
166
167            fe = FilterElement(
168                objk - config.gamma_f*thetak,
169                (1 - config.gamma_theta)*thetak)
170            filteR.addToFilter(fe)
171
172            rhok = 1 - ((theta - config.ep_i)/max(thetak, config.ep_i))
173            if rhok < config.eta1:
174                config.trust_radius = max(config.gamma_c*config.trust_radius,
175                                  config.delta_min)
176            elif rhok >= config.eta2:
177                config.trust_radius = min(config.gamma_e*config.trust_radius,
178                                  config.radius_max)
179
180            obj = problem.evaluateObj(x, y, z)
181
182            stepNorm = norm(packXYZ(x-xk, y-yk, z-zk), inf)
183            logger.setCurIter(stepNorm=stepNorm)
184
185        else:
186
187            # Solve TRSP_k
188            flag, obj = problem.TRSPk(x, y, z, xk, yk, zk,
189                                      romParam, config.trust_radius)
190            if not flag:
191                raise Exception("TRSPk fails!\n")
192
193            # Filter
194            yr = problem.evaluateDx(x)
195
196            stepNorm = norm(packXYZ(x-xk, y-yk, z-zk), inf)
197            logger.setCurIter(stepNorm=stepNorm)
198
199            theta = norm(yr - y, 1)
200            fe = FilterElement(obj, theta)
201
202            if not filteR.checkAcceptable(fe, config.theta_max) and iteration > 0:
203                logger.iterlog.rejected = True
204                config.trust_radius = max(config.gamma_c*stepNorm,
205                                  config.delta_min)
206                rebuildROM = False
207                x, y, z = cloneXYZ(xk, yk, zk)
208                continue
209
210            # Switching Condition and Trust Region update
211            if (((objk - obj) >= config.kappa_theta*
212                 pow(thetak, config.gamma_s))
213                and
214                (thetak < config.theta_min)):
215                logger.iterlog.fStep = True
216
217                config.trust_radius = min(
218                    max(config.gamma_e*stepNorm, config.trust_radius),
219                    config.radius_max)
220
221            else:
222                logger.iterlog.thetaStep = True
223
224                fe = FilterElement(
225                    obj - config.gamma_f*theta,
226                    (1 - config.gamma_theta)*theta)
227                filteR.addToFilter(fe)
228
229                # Calculate rho for theta step trust region update
230                rhok = 1 - ((theta - config.ep_i) /
231                            max(thetak, config.ep_i))
232                if rhok < config.eta1:
233                    config.trust_radius = max(config.gamma_c*stepNorm,
234                                      config.delta_min)
235                elif rhok >= config.eta2:
236                    config.trust_radius = min(
237                        max(config.gamma_e*stepNorm, config.trust_radius),
238                        config.radius_max)
239
240        # Accept step
241        rebuildROM = True
242        xk, yk, zk = cloneXYZ(x, y, z)
243        thetak = theta
244        objk = obj
245
246
247    logger.printVectors()
248#    problem.reverseTransform()
249