problem with fmin_cobyla
Hey, A bit of a newbie to this, but I have a problem which requires a dynamic set of constraints (some of which are non-linear) and I tried two versions using fmin_cobyla (both examples attached)-- one generates these constraints using lambda functions, but fmin seems to violate them (example in code). Then I tried generating them on the fly using exec on strings of functions, which observed the constraints, but failed to find the most optimal solution. (the third permutation should be higher) Can anyone help? Could not find any examples online which were more than trivial, and the docs dont seem very good. (apologies for the cryptic coding, I tried to minimize a real-world example into a shorter script) Thanx in advance, Josh ____________________________________________________________________________________ Pinpoint customers who are looking for what you sell. http://searchmarketing.yahoo.com/ import calendar import datetime import numpy from scipy.optimize import fmin_cobyla def Solve(): maxS = 17000000 dates = [datetime.date(2008,4,1),datetime.date(2008,5,1), datetime.date(2008,6,1)] vals = [.024925574678, .128886905103,.0447355248121] # this is a set of date lists, each one corresponding to one of the vals above and one of the unknown vars to solve for permutes = [[dates[1],dates[0]],[dates[2],dates[0]],[dates[2],dates[1]]] # daily limits (in the constraints, we use these at a monthly level) upLim = [114033,114033,114033] dnLim = [159646,159646,159646] allCons = genCons( dates,permutes,0,maxS,upLim,dnLim ) k = fmin_cobyla(minFunc, [5000000 for x in permutes], allCons, args=(permutes, vals),consargs=(),rhobeg=10000,rhoend=500,iprint=3,maxfun=100000) print permutes, vals # change to max and print print -1 * minFunc(k,permutes,vals) return k def minFunc( x, permutes, vals ): # we really need max, so we multiply by -1 return sum([-1*x[k]*vals[k] for k in xrange(len(permutes))]) def genCons( dates, permutes, initial, maxS, upLim, dnLim ): # generate constraints dynamically since normally we dont know how many dates there are maxSet = [] minSet = [] iSet = [] wSet = [] signSet = [] upperSet = [] lowerSet = [] fullRelList = [] fullSignList = [] for i,d in enumerate( dates ): iLimit = upLim[i] * ( calendar.mdays[d.month] + ( calendar.isleap(d.year) and d.month == 2 ) ) wLimit = dnLim[i] * ( calendar.mdays[d.month] + ( calendar.isleap(d.year) and d.month == 2 ) ) relList = [ n for n,t in enumerate( permutes ) if d in t ] signList = [ (t[0] == d and -1) or (t[1] == d and 1) for t in permutes if d in t ] #print iLimit,wLimit fullRelList.append(relList) fullSignList.append(signList) #print signList, relList iFn = lambda x: iLimit-sum([x[k]*signList[t] for t,k in enumerate(relList)]) # violated constraint!!! print iFn([6500000, 4949026, 0]) iSet.append(iFn) wFn = lambda x: sum([x[k]*signList[t] for t,k in enumerate(relList)])+wLimit wSet.append(wFn) maxFn = lambda x: maxS-initial-sum(numpy.concatenate([[x[k]*fullSignList[j][t] for t,k in enumerate(fullRelList[j])] for j in xrange(len(fullRelList))])) maxSet.append(maxFn) minFn = lambda x: sum(numpy.concatenate([[x[k]*fullSignList[j][t] for t,k in enumerate(fullRelList[j])] for j in xrange(len(fullRelList))]))-(0-initial) minSet.append(minFn) signPairs = numpy.concatenate( [ [ (x,y,relList.index(x),relList.index(y)) for y in relList if y!=x ] for x in relList ] ).tolist() for j,pair in enumerate( signPairs ): signFn = lambda x: x[pair[0]]*signList[pair[2]]*x[pair[1]]*signList[pair[3]] signSet.append(signFn) for j,item in enumerate( permutes ): iLimit = upLim[dates.index(item[0])] * (calendar.mdays[item[0].month]+(calendar.isleap(item[0].year) and item[0].month==2)) wLimit = dnLim[dates.index(item[1])] * (calendar.mdays[item[1].month]+(calendar.isleap(item[1].year) and item[1].month==2)) upperFn = lambda x: max(iLimit,wLimit) - x[j] lowerFn = lambda x: x[j] upperSet.append(upperFn) lowerSet.append(lowerFn) fullSet = iSet+wSet+maxSet+minSet+signSet+upperSet+lowerSet #print fullSet, len(fullSet) return fullSet import calendar import datetime import numpy from scipy.optimize import fmin_cobyla def Solve(): maxS = 17000000 dates = [datetime.date(2008,4,1),datetime.date(2008,5,1), datetime.date(2008,6,1)] vals = [.024925574678, .128886905103,.0447355248121] # this is a set of date lists, each one corresponding to one of the vals above and one of the unknown vars to solve for permutes = [[dates[1],dates[0]],[dates[2],dates[0]],[dates[2],dates[1]]] # daily limits (in the constraints, we use these at a monthly level) upLim = [114033,114033,114033] dnLim = [159646,159646,159646] allCons = genCons( dates,permutes,0,maxS,upLim,dnLim ) k = fmin_cobyla(minFunc, [5000000 for x in permutes], allCons, args=(permutes, vals),consargs=(),rhobeg=10000,rhoend=500,iprint=3,maxfun=100000) print permutes, vals # change to max and print print -1 * minFunc(k,permutes,vals) return k def minFunc( x, permutes, vals ): # we really need max, so we multiply by -1 return sum([-1*x[k]*vals[k] for k in xrange(len(permutes))]) def genCons( dates, permutes, initial, maxS, upLim, dnLim ): # generate constraints dynamically since normally we dont know how many dates there are maxSet = [] minSet = [] iSet = [] wSet = [] signSet = [] upperSet = [] lowerSet = [] fullRelList = [] fullSignList = [] for i,d in enumerate( dates ): iLimit = upLim[i] * ( calendar.mdays[d.month] + ( calendar.isleap(d.year) and d.month == 2 ) ) wLimit = dnLim[i] * ( calendar.mdays[d.month] + ( calendar.isleap(d.year) and d.month == 2 ) ) relList = [ n for n,t in enumerate( permutes ) if d in t ] signList = [ (t[0] == d and -1) or (t[1] == d and 1) for t in permutes if d in t ] #print iLimit,wLimit fullRelList.append(relList) fullSignList.append(signList) #print signList, relList iFn = lambda x: iLimit-sum([x[k]*signList[t] for t,k in enumerate(relList)]) # violated constraint!!! print iFn([6500000, 4949026, 0]) iSet.append(iFn) wFn = lambda x: sum([x[k]*signList[t] for t,k in enumerate(relList)])+wLimit wSet.append(wFn) maxFn = lambda x: maxS-initial-sum(numpy.concatenate([[x[k]*fullSignList[j][t] for t,k in enumerate(fullRelList[j])] for j in xrange(len(fullRelList))])) maxSet.append(maxFn) minFn = lambda x: sum(numpy.concatenate([[x[k]*fullSignList[j][t] for t,k in enumerate(fullRelList[j])] for j in xrange(len(fullRelList))]))-(0-initial) minSet.append(minFn) signPairs = numpy.concatenate( [ [ (x,y,relList.index(x),relList.index(y)) for y in relList if y!=x ] for x in relList ] ).tolist() for j,pair in enumerate( signPairs ): signFn = lambda x: x[pair[0]]*signList[pair[2]]*x[pair[1]]*signList[pair[3]] signSet.append(signFn) for j,item in enumerate( permutes ): iLimit = upLim[dates.index(item[0])] * (calendar.mdays[item[0].month]+(calendar.isleap(item[0].year) and item[0].month==2)) wLimit = dnLim[dates.index(item[1])] * (calendar.mdays[item[1].month]+(calendar.isleap(item[1].year) and item[1].month==2)) upperFn = lambda x: max(iLimit,wLimit) - x[j] lowerFn = lambda x: x[j] upperSet.append(upperFn) lowerSet.append(lowerFn) fullSet = iSet+wSet+maxSet+minSet+signSet+upperSet+lowerSet #print fullSet, len(fullSet) return fullSet def main(): return Solve()
On 07/09/07, Josh Gottlieb <yosh_6@yahoo.com> wrote:
Hey, A bit of a newbie to this, but I have a problem which requires a dynamic set of constraints (some of which are non-linear) and I tried two versions using fmin_cobyla (both examples attached)-- one generates these constraints using lambda functions, but fmin seems to violate them (example in code). Then I tried generating them on the fly using exec on strings of functions, which observed the constraints, but failed to find the most optimal solution. (the third permutation should be higher) Can anyone help? Could not find any examples online which were more than trivial, and the docs dont seem very good. (apologies for the cryptic coding, I tried to minimize a real-world example into a shorter script)
Whipping up functions on the fly, with lambda or by using def() inside a function, is a perfectly reasonable way to implement constraints. You can also use a single constraint function that takes an extra argument to tell it which constraints it should use, and pass that extra argument in through fmin_cobyla. Be warned that it's a non-trivial technical problem to implement a constrained minimizer that never evaluates the function at a point that violates the constraints; fmin_cobyla may do this, and I think the other choices in scipy may as well. Finally, remember that all the fmin_* functions are only *local* minimizers - they are supposed to find points that are local minima of the function, but if the function is not concave up, it may have many minima, and the solver doesn't even try to arrange you wind up at the lowest. You can improve your chances by starting close to the minimum, but if there are several unknown minima, you need to look for a global optimizer, which is a much more difficult problem. Scipy has a couple of rudimentary ones. Good luck, Anne M. Archibald
So what's wrong with your example? According to my results obtained from the one (1st py-file), max constraint violation is 2.37501074363e-10 Don't you forget that fmin_cobyla uses c(x)>=0 constraints? (As for me it's one more reason to use universal frameworks that cut down such edges of solvers, no needs to study each one deeply). Regards, D. Josh Gottlieb wrote:
Hey, A bit of a newbie to this, but I have a problem which requires a dynamic set of constraints (some of which are non-linear) and I tried two versions using fmin_cobyla (both examples attached)-- one generates these constraints using lambda functions, but fmin seems to violate them (example in code). Then I tried generating them on the fly using exec on strings of functions, which observed the constraints, but failed to find the most optimal solution. (the third permutation should be higher) Can anyone help? Could not find any examples online which were more than trivial, and the docs dont seem very good. (apologies for the cryptic coding, I tried to minimize a real-world example into a shorter script)
Thanx in advance, Josh
____________________________________________________________________________________ Pinpoint customers who are looking for what you sell. http://searchmarketing.yahoo.com/ ------------------------------------------------------------------------
import calendar import datetime import numpy from scipy.optimize import fmin_cobyla
def Solve(): maxS = 17000000 dates = [datetime.date(2008,4,1),datetime.date(2008,5,1), datetime.date(2008,6,1)] vals = [.024925574678, .128886905103,.0447355248121] # this is a set of date lists, each one corresponding to one of the vals above and one of the unknown vars to solve for permutes = [[dates[1],dates[0]],[dates[2],dates[0]],[dates[2],dates[1]]] # daily limits (in the constraints, we use these at a monthly level) upLim = [114033,114033,114033] dnLim = [159646,159646,159646]
allCons = genCons( dates,permutes,0,maxS,upLim,dnLim ) k = fmin_cobyla(minFunc, [5000000 for x in permutes], allCons, args=(permutes, vals),consargs=(),rhobeg=10000,rhoend=500,iprint=3,maxfun=100000) print permutes, vals # change to max and print print -1 * minFunc(k,permutes,vals) return k
def minFunc( x, permutes, vals ): # we really need max, so we multiply by -1 return sum([-1*x[k]*vals[k] for k in xrange(len(permutes))])
def genCons( dates, permutes, initial, maxS, upLim, dnLim ): # generate constraints dynamically since normally we dont know how many dates there are maxSet = [] minSet = [] iSet = [] wSet = [] signSet = [] upperSet = [] lowerSet = [] fullRelList = [] fullSignList = [] for i,d in enumerate( dates ): iLimit = upLim[i] * ( calendar.mdays[d.month] + ( calendar.isleap(d.year) and d.month == 2 ) ) wLimit = dnLim[i] * ( calendar.mdays[d.month] + ( calendar.isleap(d.year) and d.month == 2 ) ) relList = [ n for n,t in enumerate( permutes ) if d in t ] signList = [ (t[0] == d and -1) or (t[1] == d and 1) for t in permutes if d in t ] #print iLimit,wLimit fullRelList.append(relList) fullSignList.append(signList) #print signList, relList iFn = lambda x: iLimit-sum([x[k]*signList[t] for t,k in enumerate(relList)]) # violated constraint!!! print iFn([6500000, 4949026, 0]) iSet.append(iFn) wFn = lambda x: sum([x[k]*signList[t] for t,k in enumerate(relList)])+wLimit wSet.append(wFn) maxFn = lambda x: maxS-initial-sum(numpy.concatenate([[x[k]*fullSignList[j][t] for t,k in enumerate(fullRelList[j])] for j in xrange(len(fullRelList))])) maxSet.append(maxFn) minFn = lambda x: sum(numpy.concatenate([[x[k]*fullSignList[j][t] for t,k in enumerate(fullRelList[j])] for j in xrange(len(fullRelList))]))-(0-initial) minSet.append(minFn) signPairs = numpy.concatenate( [ [ (x,y,relList.index(x),relList.index(y)) for y in relList if y!=x ] for x in relList ] ).tolist() for j,pair in enumerate( signPairs ): signFn = lambda x: x[pair[0]]*signList[pair[2]]*x[pair[1]]*signList[pair[3]] signSet.append(signFn) for j,item in enumerate( permutes ): iLimit = upLim[dates.index(item[0])] * (calendar.mdays[item[0].month]+(calendar.isleap(item[0].year) and item[0].month==2)) wLimit = dnLim[dates.index(item[1])] * (calendar.mdays[item[1].month]+(calendar.isleap(item[1].year) and item[1].month==2)) upperFn = lambda x: max(iLimit,wLimit) - x[j] lowerFn = lambda x: x[j] upperSet.append(upperFn) lowerSet.append(lowerFn) fullSet = iSet+wSet+maxSet+minSet+signSet+upperSet+lowerSet #print fullSet, len(fullSet) return fullSet
------------------------------------------------------------------------
import calendar import datetime import numpy from scipy.optimize import fmin_cobyla
def Solve(): maxS = 17000000 dates = [datetime.date(2008,4,1),datetime.date(2008,5,1), datetime.date(2008,6,1)] vals = [.024925574678, .128886905103,.0447355248121] # this is a set of date lists, each one corresponding to one of the vals above and one of the unknown vars to solve for permutes = [[dates[1],dates[0]],[dates[2],dates[0]],[dates[2],dates[1]]] # daily limits (in the constraints, we use these at a monthly level) upLim = [114033,114033,114033] dnLim = [159646,159646,159646]
allCons = genCons( dates,permutes,0,maxS,upLim,dnLim ) k = fmin_cobyla(minFunc, [5000000 for x in permutes], allCons, args=(permutes, vals),consargs=(),rhobeg=10000,rhoend=500,iprint=3,maxfun=100000) print permutes, vals # change to max and print print -1 * minFunc(k,permutes,vals) return k
def minFunc( x, permutes, vals ): # we really need max, so we multiply by -1 return sum([-1*x[k]*vals[k] for k in xrange(len(permutes))])
def genCons( dates, permutes, initial, maxS, upLim, dnLim ): # generate constraints dynamically since normally we dont know how many dates there are maxSet = [] minSet = [] iSet = [] wSet = [] signSet = [] upperSet = [] lowerSet = [] fullRelList = [] fullSignList = [] for i,d in enumerate( dates ): iLimit = upLim[i] * ( calendar.mdays[d.month] + ( calendar.isleap(d.year) and d.month == 2 ) ) wLimit = dnLim[i] * ( calendar.mdays[d.month] + ( calendar.isleap(d.year) and d.month == 2 ) ) relList = [ n for n,t in enumerate( permutes ) if d in t ] signList = [ (t[0] == d and -1) or (t[1] == d and 1) for t in permutes if d in t ] #print iLimit,wLimit fullRelList.append(relList) fullSignList.append(signList) #print signList, relList iFn = lambda x: iLimit-sum([x[k]*signList[t] for t,k in enumerate(relList)]) # violated constraint!!! print iFn([6500000, 4949026, 0]) iSet.append(iFn) wFn = lambda x: sum([x[k]*signList[t] for t,k in enumerate(relList)])+wLimit wSet.append(wFn) maxFn = lambda x: maxS-initial-sum(numpy.concatenate([[x[k]*fullSignList[j][t] for t,k in enumerate(fullRelList[j])] for j in xrange(len(fullRelList))])) maxSet.append(maxFn) minFn = lambda x: sum(numpy.concatenate([[x[k]*fullSignList[j][t] for t,k in enumerate(fullRelList[j])] for j in xrange(len(fullRelList))]))-(0-initial) minSet.append(minFn) signPairs = numpy.concatenate( [ [ (x,y,relList.index(x),relList.index(y)) for y in relList if y!=x ] for x in relList ] ).tolist() for j,pair in enumerate( signPairs ): signFn = lambda x: x[pair[0]]*signList[pair[2]]*x[pair[1]]*signList[pair[3]] signSet.append(signFn) for j,item in enumerate( permutes ): iLimit = upLim[dates.index(item[0])] * (calendar.mdays[item[0].month]+(calendar.isleap(item[0].year) and item[0].month==2)) wLimit = dnLim[dates.index(item[1])] * (calendar.mdays[item[1].month]+(calendar.isleap(item[1].year) and item[1].month==2)) upperFn = lambda x: max(iLimit,wLimit) - x[j] lowerFn = lambda x: x[j] upperSet.append(upperFn) lowerSet.append(lowerFn) fullSet = iSet+wSet+maxSet+minSet+signSet+upperSet+lowerSet #print fullSet, len(fullSet) return fullSet
def main(): return Solve()
------------------------------------------------------------------------
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
participants (3)
-
Anne Archibald -
dmitrey -
Josh Gottlieb