[solved] Goal seek macro

Need help, or want to share a macro? Post here!
Forum rules
Be nice to others! Respect the FreeCAD code of conduct!
User avatar
Chris_G
Veteran
Posts: 2579
Joined: Tue Dec 31, 2013 4:10 pm
Location: France
Contact:

Re: [Help] Goal seek macro

Post by Chris_G »

davidosterberg wrote: Sat Jan 23, 2021 3:42 pm Scipy is awesome, definitely the right way. You can even tell it to respect bounds on the parameters (i.e max and min value)
by changing method to e.g "L-BFGS-B" and setting the bounds parameter to a list of (min, max) tuples.
Thanks.
I was interested with parameters bounds, but get lost by all available algorithms.
edwilliams16 wrote: Sat Jan 23, 2021 10:49 pm The minimization may work better using

Code: Select all

return (o2.ddGoal.Value - o2.ddResult.Value)**2
because the latter function is smooth at the minimum, and gradient descent methods will automatically take smaller steps as the minimum is approached.
Thanks, that's interesting. I'll give a try at that.
kisolre
Veteran
Posts: 4162
Joined: Wed Nov 21, 2018 1:13 pm

Re: [Help] Goal seek macro

Post by kisolre »

Chris_G wrote: Sun Jan 24, 2021 8:29 am.
edwilliams16 wrote: Sat Jan 23, 2021 10:49 pm.
davidosterberg wrote: Sat Jan 23, 2021 3:42 pm.
openBrain wrote: Sat Jan 23, 2021 12:19 am.
Sorry for the bump but can you suggest a way to exit the macro on failed Recompute? I tried generating an Exception but that just prints the message in report view and minimize happily continues until 100 iterations pass. Which takes time without doing any meaningful job.
This version (when the minimize starts from 1) causes the sketch in GoalSeekTest003 to fail and so does the Recompute:

Code: Select all

import FreeCAD
from scipy.optimize import minimize

sel = Gui.Selection.getSelectionEx()
if len(sel)!=2:
    raise Exception("Select 2 DD objects - Parameter and Goal")
o1=sel[0].Object
o2=sel[1].Object
def compute(value):
    o1.ddParameter = value[0] 
    print(o1.ddParameter)
    if not(App.ActiveDocument.recompute()):
       raise Exception("Recompuute failed, Please check model for errors (broken sketckes, toponaming, ... ")

    return (o2.ddGoal.Value - o2.ddResult.Value)**2

#res = minimize(compute, o1.ddParameter, method='Nelder-Mead', options={"maxiter": 100, "disp": True})
res = minimize(compute, 1, method='Nelder-Mead', options={"maxiter": 100, "disp": True})

print("For Parameter = {}, Result = {}".format(o1.ddParameter, o2.ddResult))
davidosterberg
Posts: 529
Joined: Fri Sep 18, 2020 5:40 pm

Re: [Help] Goal seek macro

Post by davidosterberg »

kisolre wrote: Sun Jan 24, 2021 9:44 am Sorry for the bump but can you suggest a way to exit the macro on failed Recompute?
There is no proper way (cf https://github.com/scipy/scipy/issues/7306).

I suggest

Code: Select all

return float('nan')
It should make the optimization stop, but I havent tried it.
keithsloan52
Veteran
Posts: 2756
Joined: Mon Feb 27, 2012 5:31 pm

Re: [Help] Goal seek macro

Post by keithsloan52 »

Chris_G wrote: Sat Jan 23, 2021 9:30 am If you can install python.scipy, this little script seems to work :

Code: Select all

import FreeCAD
from scipy.optimize import minimize

doc1 = FreeCAD.getDocument('GoalSeekTest001')
o1 = doc1.getObject('dd')
o2 = doc1.getObject('dd001')


def compute(value):
    o1.ddParameter = value[0]
    doc1.recompute()
    return abs(o2.ddGoal.Value - o2.ddResult.Value)

res = minimize(compute, 1, method='Nelder-Mead', options={"maxiter": 2000, "disp": True})

print("For Height = {}, Volume = {}".format(o1.ddParameter, o2.ddResult))
Minor query : Would it not br better to use o1.recompute() rather than a recompute of the whole document, which may contain a large number of objects?
openBrain
Veteran
Posts: 9034
Joined: Fri Nov 09, 2018 5:38 pm
Contact:

Re: [Help] Goal seek macro

Post by openBrain »

keithsloan52 wrote: Sun Jan 24, 2021 10:37 am Minor query : Would it not br better to use o1.recompute() rather than a recompute of the whole document, which may contain a large number of objects?
Nope. 'o1' is just the parameter. You need to at least recompute the result in 'o2'.
Maybe doing 'o2.recompute(True)' -- 'True' is for recursive because you also need intermediate objects to be recomputed -- would be enough, but I'm not sure if recursive recomputing took into account the expression bindings.
kisolre
Veteran
Posts: 4162
Joined: Wed Nov 21, 2018 1:13 pm

Re: [Help] Goal seek macro

Post by kisolre »

Code: Select all

if not(App.ActiveDocument.recompute()):
This actually does not work because Recompute returns number of objects recomputed. What I need to check is if the document is still touched after the recompute, which should be the state if recompute does not complete (fail). But I cant find how to do that.
The API docs say

Code: Select all

App::Document::isTouched()
◆ isTouched()
bool Document::isTouched 	(	void 		)	 const
check if there is any touched object in this document
References App::DocumentP::objectArray.
but trying that in the console gives

Code: Select all

>>> App.ActiveDocument.isTouched()
Traceback (most recent call last):
  File "<input>", line 1, in <module>
AttributeError: 'App.Document' object has no attribute 'isTouched'
>>> 
I surely dont understand much Python but is this wrong?
openBrain
Veteran
Posts: 9034
Joined: Fri Nov 09, 2018 5:38 pm
Contact:

Re: [Help] Goal seek macro

Post by openBrain »

kisolre wrote: Sun Jan 24, 2021 1:00 pm I surely dont understand much Python but is this wrong?
I didn't find any trace of the function you pointed out.
There is 'DocumentObject::isTouched()' but as it says, it applies to objects, not directly to document.

You can use an iterator to test if there is still touched objects in the document (snippets returns 'True' if at least one object is touched, 'False' otherwise):

Code: Select all

next((True for o in App.ActiveDocument.Objects if o.State==['Touched']),False)
Another solution would be to recompute a second time the document. If output isn't 0, there is an error in the document.

HTH
kisolre
Veteran
Posts: 4162
Joined: Wed Nov 21, 2018 1:13 pm

Re: [Help] Goal seek macro

Post by kisolre »

openBrain wrote: Sun Jan 24, 2021 1:19 pm I didn't find any trace of the function you pointed out.
Here is a link https://freecad.github.io/SourceDoc/d8/ ... 6bdfca2af0
openBrain wrote: Sun Jan 24, 2021 1:19 pm Another solution would be to recompute a second time the document. If output isn't 0, there is an error in the document.
This actually works. and then the exception properly kills the minimize and the macro.

So here is the latest "complete" version:

Code: Select all

# Goal Seeek FreeCAD macro
# Author: Chris_G ,Kisolre
# License: LGPL-3.0-or-later
#
# Usage: This macro uses scipy and DynamicData WB to adjust a Parameter so that Result reaches Goal.
#
# An example is adjusting the height of water glass model so the internal volume reaches desired value (260ml for example)
#
# For entering desired values two DD objects are used. 
# First one should have a property named "ddParameter" under it which should be used by expression in the model (the glass height in the above example).
# Second one should have properties named "ddResult" and "ddGoal" under it.
# ddResult should reference through expressions some value from the model dependent on ddParameter (volume of the glass (Body.Shape.Volume)).
# ddGoal is the desired value that ddResult should reach (260ml)
# 
# Ctrl-select the dd objects in proper order, execute macro, wait for result. Each iteration recomputes the documend so depending on the model this can be slow.

import FreeCAD
from scipy.optimize import minimize

sel = Gui.Selection.getSelectionEx()
if len(sel)!=2:
    raise Exception("Select 2 DD objects - Parameter and Goal")
o1=sel[0].Object
o2=sel[1].Object
def compute(value):
    o1.ddParameter = value[0] 
    print(o1.ddParameter)
    App.ActiveDocument.recompute()
    # This will check if there are still objects to recompute - previous recompute failed.
    if App.ActiveDocument.recompute()!=0: 
        raise Exception ("Recompuute failed, Please check model for errors (broken sketckes, toponaming, ... ")

    return (o2.ddGoal.Value - o2.ddResult.Value)**2

res = minimize(compute, o1.ddParameter, method='Nelder-Mead', options={"maxiter": 100, "disp": True})

print("For Parameter = {}, Result = {}".format(o1.ddParameter, o2.ddResult))
openBrain
Veteran
Posts: 9034
Joined: Fri Nov 09, 2018 5:38 pm
Contact:

Re: [Help] Goal seek macro

Post by openBrain »

kisolre wrote: Sun Jan 24, 2021 2:08 pm
openBrain wrote: Sun Jan 24, 2021 1:19 pm I didn't find any trace of the function you pointed out.
Here is a link https://freecad.github.io/SourceDoc/d8/ ... 6bdfca2af0
OK. I missed the point. This hasn't been exposed to Python API ATM, so available only in C++.
openBrain
Veteran
Posts: 9034
Joined: Fri Nov 09, 2018 5:38 pm
Contact:

Re: [Help] Goal seek macro

Post by openBrain »

As a hobby, below a quick and dirty seeking function using basic dichotomy

Code: Select all

SEEKMAX = 100 # max number of iterations to found a range where a solution exists
DICHOMAX = 100 # max number of iterations to found a solution inside the found range
SEEKFACTOR = 0.5 # range size around nominal value, shall be in ]0, 1[ ; 0.5 typical ; reducing may improve finding a solution on periodic functions
DICHOFACTOR = 0.5 # dichotomy factor, no valuable reason to set something else than 0.5
THRESHOLD = 1e-7 # deviation with goal at which a solution is considered found

from PySide2 import QtWidgets
import FreeCAD as App
import FreeCADGui as Gui

doc = App.ActiveDocument
ddp = doc.dd
ddg = doc.dd001
goal = ddg.ddGoal

sb = Gui.getMainWindow().findChild(QtWidgets.QProgressBar)

def recompute(val = ddp.ddParameter):
    if not ddg.recompute(True):
        raise Exception("Recompute error")

def goalseek(baseValue, currValue = None, itseek = 1, itdicho = 1, up = False, down = False):
    global SEEKFACTOR
    sb.setValue((itseek+itdicho)/(SEEKMAX+DICHOMAX)*100)
    Gui.updateGui()
    if itseek > SEEKMAX+1:
        App.Console.PrintWarning("Maximum number of seek iterations reached without solution range\n")
        return None
    if itdicho > DICHOMAX+1:
        App.Console.PrintWarning("Maximum number of dichotomy iterations reached without solution\n")
        return None
    if currValue:
            ddp.ddParameter = currValue
    else:
            ddp.ddParameter = baseValue
    recompute() 
    currResult = ddg.ddResult
    ## SEEK phase
    if up == down:
        upValue = baseValue + baseValue*SEEKFACTOR
        ddp.ddParameter = upValue
        recompute()
        upResult = ddg.ddResult
        downValue = baseValue - baseValue*SEEKFACTOR
        ddp.ddParameter = downValue
        recompute()
        downResult = ddg.ddResult
        if (currResult < goal) == (upResult > goal): # found a range where a solution exists upside
            return goalseek(baseValue, baseValue, itseek, itdicho, True, False)
        if (currResult > goal) == (downResult < goal): # found a range where a solution exists downside
            return goalseek(baseValue, baseValue, itseek, itdicho, False, True)
        if (upResult > currResult) == (currResult > downResult) == (downResult > goal): # found a linear range, solution shoud be downside
            return goalseek(downValue, currValue, itseek+1, itdicho, up, down)
        if (upResult > currResult) == (currResult > downResult) == (upResult < goal): # found a linear range, solution shoud be upside
            return goalseek(upValue, currValue, itseek+1, itdicho, up, down)
        ## found a concave/convexe range
        upDiff = upResult - currResult
        downDiff = currResult - downResult
        SEEKFACTOR = SEEKFACTOR / 2
        if ((currResult < goal) == (downDiff < 0)) != (abs(downDiff) > abs(upDiff)):
            return goalseek(baseValue, currValue, itseek+1, itdicho, up, down)
        else:
            return goalseek(baseValue, currValue, itseek+1, itdicho, up, down)
    ## DICHO phase
    if up:
        currValue = currValue + baseValue*SEEKFACTOR*DICHOFACTOR**itdicho
    if down:
        currValue = currValue - baseValue*SEEKFACTOR*DICHOFACTOR**itdicho
    ddp.ddParameter = currValue
    recompute()
    if abs(ddg.ddResult - goal) <= THRESHOLD:
        print(f"""Found a solution with value : {currValue}
            Deviation : {abs(ddg.ddResult - goal)}
            Seek iterations : {itseek}
            Dicho iterations : {itdicho}""")
        return currValue
    if (currResult < goal) == (ddg.ddResult < goal):
        return goalseek(baseValue, currValue, itseek, itdicho+1, up, down)
    else:
        return goalseek(baseValue, currValue, itseek, itdicho+1, not up, not down)

def runSeek():
    App.Console.PrintMessage(f"Seeking for goal : {goal}\n")
    sb.setValue(0)
    sb.show()
    startVal = ddp.ddParameter
    res = None
    try:
        res = goalseek(startVal)
    except Exception:
        App.Console.PrintError(f"Error : can't recompute with value : {ddp.ddParameter}\n")
    if not res:
        ddp.ddParameter = startVal # Restore start value if no solution found
    sb.hide()

runSeek()
I didn't test it extensively but it should deal quite well with any polynomial (and maybe even periodic) transfer function. :)
Post Reply