L. Jaulin (2021). A boundary approach for set inversion, Engineering Applications of Artificial Intelligence, Volume 100, April 2021

Code

from pyibex import *
from pyibex.geometry import SepPolygon
from vibes import vibes
from scipy.optimize import linprog

dd = Function('dd1', 'dd2', 'max((dd1-8)^2-0.001^2,(dd2-4)^2-0.001^2)')
Cdd = CtcFwdBwd(dd)

fdy= Function('d1','d2','d3','y1','y2','(d2-d1-y1;d3-d2-y2;d3-d1-y1-y2)')
Cdy=CtcFwdBwd(fdy)

fxd=Function('x1', 'x2','d1','d2','d3',
        '(sqrt((4-x1)^2+(6-x2)^2)-d1;sqrt((13-x1)^2+(7-x2)^2)-d2;sqrt((16-x1)^2+(10-x2)^2)-d3)')
Cxd=CtcFwdBwd(fxd)


class myCtc(Ctc):
    def __init__(C):
        Ctc.__init__(C, 2)
    def contract(C, X):
        x1, x2 = X[0], X[1]

        d1=Interval(-oo,oo)
        d2=Interval(-oo,oo)
        d3=Interval(-oo,oo)
        y1=Interval(-oo,oo)
        y2=Interval(-oo,oo)

        for i in range(1,10):
            XD=IntervalVector([x1,x2,d1,d2,d3])
            Cxd.contract(XD)
            x1,x2,d1,d2,d3=XD[0],XD[1],XD[2],XD[3],XD[4]

            DY=IntervalVector([d1,d2,d3,y1,y2])
            Cdy.contract(DY)
            d1, d2, d3, y1, y2=DY[0],DY[1],DY[2],DY[3],DY[4]

            DD=IntervalVector([y1,y2])
            Cdd.contract(DD)
            y1, y2=DD[0],DD[1]

            #print("d1, d2, d3",d1, d2, d3)
            #print("y1, y2",y1, y2)

            #Da = IntervalVector([d1, d2, d3])
            #Ya = IntervalVector([y1, y2])
            #Db,Yb=contract_polytop(Da, Ya)  Dans notre cas, ne contracte rien du tout.
            #d1, d2, d3=Db[0],Db[1],Db[2]
            #print('Db=',Db)

            DY=IntervalVector([d1,d2,d3,y1,y2])
            Cdy.contract(DY)
            d1, d2, d3, y1, y2=DY[0],DY[1],DY[2],DY[3],DY[4]
            #print(d1,d2,d3)

            XD=IntervalVector([x1,x2,d1,d2,d3])
            Cxd.contract(XD)
            x1,x2,d1,d2,d3=XD[0],XD[1],XD[2],XD[3],XD[4]

        X[0]=x1
        X[1]=x2

C2= myCtc()

def contract_polytop(X,Y):
    A = [[-1,1,0,-1,0],[0,-1,1,0,-1]]
    b = [0, 0]
    bb=[(X[0].lb(), X[0].ub()), (X[1].lb(), X[1].ub()),  (X[2].lb(), X[2].ub()), (Y[0].lb(), Y[0].ub()), (Y[1].lb(), Y[1].ub())]
    n=5
    Z=[[0,0],[0,0],[0,0],[0,0],[0,0]]
    for i in range(0,n):
        p=[0.,0.]
        for j in [0,1]:
            c = [0] * n
            c[i]=(j-0.5)*2.0
            res=linprog(c, A_eq=A, b_eq=b, bounds=bb)
            p[j]=res["fun"]
        Z[i]=Interval(p[1],-p[0])
    X=[Interval(Z[0]),Interval(Z[1]),Interval(Z[2])]
    Y=[Interval(Z[3]),Interval(Z[4])]
    return X,Y


X0 = IntervalVector([[0.94,1.06],[1.94,2.06]])

vibes.beginDrawing()
color = {'color_in':'black[red]', 'color_out':'blue[cyan]', 'color_maybe':'white[yellow]'}
color = {'color_in':'gray[white]', 'color_out':'gray[white]', 'color_maybe':'black[yellow]'}
pySIVIA(X0, C2, 0.0001, **color)
vibes.axisAuto()
vibes.setFigureSize(500,500)
vibes.endDrawing()