|
from __future__ import print_function
from __future__ import absolute_import
from __future__ import division
from builtins import range
from builtins import object
from gpkit import Model, Variable, SignomialsEnabled
from gpkit.nomials import SignomialInequality, MonomialEquality
from gpkit.exceptions import InvalidGPConstraint
from gpkit.nomials import SingleSignomialEquality
import numpy as np
from time import time
import warnings
from scipy.stats import norm
from .robust_gp_tools import RobustGPTools
from .equivalent_posynomials import EquivalentPosynomials
from .equivalent_models import TwoTermBoydModel
from .twoterm_approximation import TwoTermApproximation
from .robustify_large_posynomial import RobustifyLargePosynomial
from .linearize_twoterm_posynomials import LinearizeTwoTermPosynomials
class RobustnessSetting(object):
def __init__(self, **options):
self._options = {
'gamma': 1,
'simpleModel': False,
'numberOfRegressionPoints': 2,
'numberOfRegressionPointsElliptical': 25,
'linearizeTwoTerm': True,
'enableSP': True,
'boyd': False,
'twoTerm': True,
'simpleTwoTerm': False,
'smartTwoTermChoose': False,
'allowedNumOfPerms': 30,
'linearizationTolerance': 0.001,
'minNumOfLinearSections': 12,
'maxNumOfLinearSections': 20,
'iterationsRelativeTolerance': 1e-4,
'iterationLimit': 10,
'probabilityOfSuccess': 0.9,
'lognormal': True
}
for key, value in options.items():
self._options[key] = value
if self._options['twoTerm']:
self._options['linearizeTwoTerm'] = True
self._options['enableSP'] = False
if self._options['simpleModel']:
self._options['allowedNumOfPerms'] = 1
def get(self, option_name):
return self._options[option_name]
def set(self, option_name, value):
self._options[option_name] = value
class RobustModel(object):
"""
RobustModel extends gpkit.Model through the robust counterpart.
It uses the nominal solution of the GP or SP to
"""
def __init__(self, model, type_of_uncertainty_set, **options):
self.nominal_model = model
self.substitutions = model.substitutions
self.type_of_uncertainty_set = type_of_uncertainty_set
self.setting = RobustnessSetting(**options)
slopes_intercepts = LinearizeTwoTermPosynomials. \
two_term_posynomial_linearization_coeff(self.setting.get('minNumOfLinearSections'))
self.robust_solve_properties = {'setuptime': 0,
'numoflinearsections': self.setting.get('minNumOfLinearSections'),
'slopes': slopes_intercepts[0],
'intercepts': slopes_intercepts[1]
}
self.number_of_stds = norm.ppf(self.setting.get("probabilityOfSuccess") / 2.0 + 0.5)
if 'nominalsolve' in options:
self.nominal_solve = options['nominalsolve']
else:
self.nominal_solve = RobustModel.internalsolve(model, verbosity=0)
self.nominal_solution = self.nominal_solve.get('variables')
self.nominal_cost = self.nominal_solve['cost']
self._sequence_of_rgps = []
self._robust_model = None
self.lower_approximation_is_feasible = False
if self.type_of_uncertainty_set == 'box':
self.dependent_uncertainty_set = False
else:
self.dependent_uncertainty_set = True
if self.type_of_uncertainty_set == 'elliptical':
self.setting.set('numberOfRegressionPoints', self.setting.get('numberOfRegressionPointsElliptical'))
self.ready_gp_constraints = []
self.to_linearize_gp_posynomials = []
self.large_gp_posynomials = []
self.sp_constraints = []
self.sp_equality_constraints = []
equality_constraints = False
if self.setting.get('boyd'):
self.setting.set('iterationLimit', 1)
try:
safe_model = TwoTermBoydModel(model)
except InvalidGPConstraint:
raise Exception("boyd's formulation is not supported for sp models")
safe_model_constraints = safe_model.flat(constraintsets=False)
del safe_model
for cs in safe_model_constraints:
if isinstance(cs, MonomialEquality):
self.ready_gp_constraints += [cs]
equality_constraints = True
else:
p = cs.as_posyslt1()[0]
if len(p.exps) == 1:
robust_monomial = self.robustify_monomial(p)
self.ready_gp_constraints += [robust_monomial <= 1]
else:
self.to_linearize_gp_posynomials += [p]
del safe_model_constraints
if equality_constraints:
warnings.warn('equality constraints will not be robustified')
self.number_of_gp_posynomials = 0
return
all_constraints = model.flat(constraintsets=False)
gp_posynomials = []
for cs in all_constraints:
if isinstance(cs, SingleSignomialEquality):
self.sp_equality_constraints.append(cs)
elif isinstance(cs, SignomialInequality):
self.sp_constraints.append(cs)
elif isinstance(cs, MonomialEquality):
self.ready_gp_constraints += [cs]
equality_constraints = True
else:
gp_posynomials += cs.as_posyslt1()
self.number_of_gp_posynomials = len(gp_posynomials)
constraints_posynomials_tuple = self.classify_gp_constraints(gp_posynomials)
self.ready_gp_constraints += constraints_posynomials_tuple[0]
self.to_linearize_gp_posynomials = constraints_posynomials_tuple[1]
self.large_gp_posynomials = constraints_posynomials_tuple[2]
if equality_constraints:
warnings.warn('equality constraints will not be robustified')
def setup(self, verbosity=0, **options):
for option, key in options.items():
self.setting.set(option, key)
start_time = time()
old_solution = self.nominal_solve
reached_feasibility = 0
for count in range(self.setting.get('iterationLimit')):
if verbosity > 0:
print("iteration %s" % (count + 1))
ready_sp_constraints, to_linearize_sp_posynomials, large_sp_posynomials = self. \
approximate_and_classify_sp_constraints(old_solution, self.number_of_gp_posynomials)
ready_constraints = self.ready_gp_constraints + ready_sp_constraints
to_linearize_posynomials = self.to_linearize_gp_posynomials + to_linearize_sp_posynomials
large_posynomials = self.large_gp_posynomials + large_sp_posynomials
permutation_indices = self.new_permutation_indices(old_solution, large_posynomials)
two_term_data_posynomials = []
for i, two_term_approximation in enumerate(large_posynomials):
permutation = two_term_approximation.list_of_permutations[permutation_indices[i]]
no_data, data = TwoTermApproximation. \
two_term_equivalent_posynomial(two_term_approximation.p, i, permutation, False)
ready_constraints += no_data
two_term_data_posynomials += [constraint.as_posyslt1()[0] for constraint in data]
two_term_data_posynomials += to_linearize_posynomials
if reached_feasibility:
self._robust_model, _ = self. \
linearize_and_return_upper_lower_models(two_term_data_posynomials,
self.robust_solve_properties['numoflinearsections'],
ready_constraints)
new_solution = RobustModel.internalsolve(self._robust_model, verbosity=0)
else:
try:
self.robust_solve_properties['numoflinearsections'], new_solution, self._robust_model = self. \
find_number_of_piece_wise_linearization(two_term_data_posynomials, ready_constraints)
reached_feasibility += 1
except Exception:
self.robust_solve_properties['numoflinearsections'], new_solution, self._robust_model = self. \
find_number_of_piece_wise_linearization(two_term_data_posynomials, ready_constraints,
feasible=True)
rel_tol = np.abs((new_solution['cost'] - old_solution['cost']) / old_solution['cost'])
if verbosity > 0:
if not reached_feasibility:
print("feasibility is not reached yet")
elif reached_feasibility == 1:
print("feasibility is reached")
print("relative tolerance = %s" % rel_tol)
if reached_feasibility <= 1 and two_term_data_posynomials:
self.robust_solve_properties['slopes'], self.robust_solve_properties['intercepts'], _, _, _ = \
LinearizeTwoTermPosynomials.two_term_posynomial_linearization_coeff(
self.robust_solve_properties['numoflinearsections'])
self._sequence_of_rgps.append(self._robust_model)
if rel_tol <= self.setting.get('iterationsRelativeTolerance'):
break
else:
old_solution = new_solution
if reached_feasibility < 1:
raise Exception("feasibility is not reached. If the solution seems to converge, try "
"increasing iterationLimit = %s. Increasing the allowed number of permutations might also "
"help" % self.setting.get('iterationLimit'))
self.robust_solve_properties['setuptime'] = time() - start_time
def robustsolve(self, verbosity=1, **options):
if self._robust_model is None:
self.setup(verbosity, **options)
try:
sol = self._robust_model.solve(verbosity=verbosity)
except InvalidGPConstraint:
sol = self._robust_model.localsolve(verbosity=verbosity)
if verbosity > 0:
print ("solving needed %s iterations." % len(self._sequence_of_rgps))
print ("setting up took %s seconds." % self.robust_solve_properties['setuptime'])
sol.update(self.robust_solve_properties)
return sol
def approximate_and_classify_sp_constraints(self, solution, number_of_gp_posynomials):
sp_gp_approximation = []
with SignomialsEnabled():
for cs in self.sp_constraints:
css = SignomialInequality(cs.left.sub(solution["constants"]), cs.oper, cs.right.sub(solution["constants"]))
sp_gp_approximation.append(css.as_gpconstr(x0=solution["freevariables"]).as_posyslt1()[0])
ready_sp_constraints, to_linearize_sp_posynomials, large_sp_posynomial = self.\
classify_gp_constraints(sp_gp_approximation, number_of_gp_posynomials)
for cs in self.sp_equality_constraints:
css = SingleSignomialEquality(cs.left.sub(solution["constants"]), cs.right.sub(solution["constants"]))
ready_sp_constraints.append(css.as_gpconstr(x0=solution["freevariables"]))
return ready_sp_constraints, to_linearize_sp_posynomials, large_sp_posynomial
def classify_gp_constraints(self, gp_posynomials, offset=0):
data_gp_posynomials = []
ready_gp_constraints = []
for i, p in enumerate(gp_posynomials):
equivalent_p = EquivalentPosynomials(p, i + offset, self.setting.get('simpleModel'),
self.dependent_uncertainty_set)
no_data, data = equivalent_p.no_data_constraints, equivalent_p.data_constraints
data_gp_posynomials += [posy.as_posyslt1()[0] for posy in data]
ready_gp_constraints += no_data
to_linearize_gp_posynomials = []
large_gp_posynomials = []
for i, p in enumerate(data_gp_posynomials):
if len(p.exps) == 1:
robust_monomial = self.robustify_monomial(p)
ready_gp_constraints += [robust_monomial <= 1]
elif len(p.exps) == 2 and self.setting.get('linearizeTwoTerm'):
to_linearize_gp_posynomials += [p]
else:
if self.setting.get('twoTerm'):
two_term_approximation = TwoTermApproximation(p, self.setting)
large_gp_posynomials.append(two_term_approximation)
else:
robust_large_p = RobustifyLargePosynomial(p, self.type_of_uncertainty_set,
self.number_of_stds, self.setting)
ready_gp_constraints += robust_large_p. \
robustify_large_posynomial(self.type_of_uncertainty_set, i + offset, self.setting)
return ready_gp_constraints, to_linearize_gp_posynomials, large_gp_posynomials
def robustify_monomial(self, monomial):
new_monomial_exps = RobustGPTools. \
only_uncertain_vars_monomial(monomial.exps[0])
m_direct_uncertain_vars = [var for var in list(new_monomial_exps.keys()) if RobustGPTools.is_uncertain(var)]
l_norm = 0
for var in m_direct_uncertain_vars:
eta_min, eta_max = RobustGPTools.generate_etas(var)
scale = eta_max
exponent = -new_monomial_exps.get(var.key)
pert = exponent * scale
if self.type_of_uncertainty_set == 'box':
l_norm += np.abs(pert)
elif self.type_of_uncertainty_set == 'elliptical':
l_norm += pert ** 2
elif self.type_of_uncertainty_set == 'one norm':
l_norm = max(l_norm, np.abs(pert))
else:
raise Exception('This type of set is not supported')
if self.type_of_uncertainty_set == 'elliptical':
l_norm = np.sqrt(l_norm)
g = self.setting.get('gamma')
# Fifth order Taylor approx of the e**gamma, so that gamma can be a variable
robust_monomial = monomial * (1.+g+1./2.*g**2+1./6.*g**3+1./24.*g**4+1./120.*g**5)**l_norm
return robust_monomial
def robustify_set_of_monomials(self, set_of_monomials, feasible=False):
robust_set_of_monomial_constraints = []
slackvar = Variable()
for monomial in set_of_monomials:
robust_monomial = self.robustify_monomial(monomial)
robust_set_of_monomial_constraints += [robust_monomial <= slackvar ** feasible]
robust_set_of_monomial_constraints += [slackvar >= 1, slackvar <= 1000]
return robust_set_of_monomial_constraints, slackvar
def calculate_value_of_two_term_approximated_posynomial(self, two_term_approximation, index_of_permutation,
solution):
permutation = two_term_approximation.list_of_permutations[index_of_permutation]
number_of_two_terms = int(len(permutation) / 2)
num_of_linear_sections = self.robust_solve_properties['numoflinearsections']
slopes = self.robust_solve_properties['slopes']
intercepts = self.robust_solve_properties['intercepts']
values = []
mons = two_term_approximation.chop()
for i in range(number_of_two_terms):
monomials = []
first_monomial = mons[2*i]
second_monomial = mons[2*i+1]
monomials += [first_monomial]
for j in range(num_of_linear_sections - 2):
monomials += [first_monomial ** slopes[num_of_linear_sections - 3 - j] *
second_monomial ** slopes[j] * np.exp(intercepts[j])]
monomials += [second_monomial]
subs_monomials = []
for j in range(len(monomials)):
# st3 = time()
robust_monomial = self.robustify_monomial(monomials[j])
monomials[j] = robust_monomial.sub(solution['variables'])
# print "subs for a monomial is taking too much time", time()-st3
subs_monomials.append(monomials[j].cs[0])
values.append(max(subs_monomials))
if number_of_two_terms % 2 != 0:
monomial = mons[len(permutation) - 1]
robust_monomial = self.robustify_monomial(monomial)
monomial = robust_monomial.sub(solution['variables'])
values.append(monomial.cs[0])
return sum(values)
def find_permutation_with_minimum_value(self, two_term_approximation, solution):
minimum_value = np.inf
minimum_index = len(two_term_approximation.list_of_permutations)
for i in range(len(two_term_approximation.list_of_permutations)):
temp_value = self. \
calculate_value_of_two_term_approximated_posynomial(two_term_approximation, i, solution)
if temp_value < minimum_value:
minimum_value = temp_value
minimum_index = i
return minimum_index
def linearize_and_return_upper_lower_models(self, two_term_data_posynomials, r, ready_constraints, feasible=False):
no_data_upper_constraints = []
no_data_lower_constraints = []
data_posynomials = []
for i, two_term_p in enumerate(two_term_data_posynomials):
linearize_p = LinearizeTwoTermPosynomials(two_term_p)
no_data_upper, no_data_lower, data = linearize_p. \
linearize_two_term_posynomial(i, r)
no_data_upper_constraints += no_data_upper
no_data_lower_constraints += no_data_lower
data_posynomials += [constraint.as_posyslt1()[0] for constraint in data]
del linearize_p, no_data_lower, no_data_upper
data_constraints, slackvar = self.robustify_set_of_monomials(data_posynomials, feasible)
upper_cons, lower_cons = [no_data_upper_constraints, ready_constraints, data_constraints], \
[no_data_lower_constraints, ready_constraints, data_constraints]
model_upper = Model(self.nominal_model.cost * slackvar ** (100 * feasible), upper_cons)
model_lower = Model(self.nominal_model.cost * slackvar ** (100 * feasible), lower_cons)
model_upper.substitutions.update(self.substitutions)
model_lower.substitutions.update(self.substitutions)
model_upper.unique_varkeys, model_lower.unique_varkeys = [self.nominal_model.varkeys] * 2
model_upper.reset_varkeys()
model_lower.reset_varkeys()
del upper_cons, lower_cons, no_data_lower_constraints, no_data_upper_constraints, data_posynomials
return model_upper, model_lower
def find_number_of_piece_wise_linearization(self, two_term_data_posynomials, ready_constraints, feasible=False):
the_min_r = self.setting.get('minNumOfLinearSections')
the_max_r = self.setting.get('maxNumOfLinearSections')
r = None
error = None
sol_upper = None
model_upper = None
while the_min_r <= the_max_r:
r = int((the_min_r + the_max_r) / 2.0)
model_upper, model_lower = self. \
linearize_and_return_upper_lower_models(two_term_data_posynomials, r, ready_constraints, feasible)
upper_model_infeasible = 0
try:
sol_upper = RobustModel.internalsolve(model_upper, verbosity=0)
except RuntimeWarning:
upper_model_infeasible = 1
try:
sol_lower = RobustModel.internalsolve(model_lower, verbosity=0)
except RuntimeWarning:
raise Exception("The model is infeasible")
if not two_term_data_posynomials:
self.robust_solve_properties['upperLowerRelError'] = 0
return 0, sol_upper, model_upper
if upper_model_infeasible != 1:
error = (sol_upper.get('cost') - sol_lower.get('cost')) / sol_lower.get('cost')
if error <= self.setting.get('linearizationTolerance'):
the_max_r = r
else:
the_min_r = r + 1
elif r == self.setting.get('maxNumOfLinearSections'):
self.lower_approximation_is_feasible = True
raise Exception("The model is infeasible. The lower approximation of the model is feasible, try "
"increasing the maximum number of linear sections")
else:
the_min_r = r + 1
del model_lower, sol_lower
if the_max_r == the_min_r and r == the_max_r:
break
self.robust_solve_properties['upperLowerRelError'] = error
return r, sol_upper, model_upper
def new_permutation_indices(self, solution, large_posynomials):
permutation_indices = []
for two_term_approximation in large_posynomials:
permutation_indices.append(self.find_permutation_with_minimum_value(two_term_approximation, solution))
return permutation_indices
@staticmethod
def internalsolve(model, verbosity=0):
try:
return model.solve(verbosity=verbosity)
except InvalidGPConstraint:
return model.localsolve(verbosity=verbosity)
def get_robust_model(self):
if self.sp_constraints:
return self._sequence_of_rgps
else:
return self._robust_model
def nominalsolve(self):
return self.nominal_solve
|
