Optimization in pyomo.dae

[Jinhyeun Kim]

Hello.

I have a question regarding the optimization in pyomo.dae.

I discretized the equation before using the m.obj = Objective(expr=1) and saw how it works.

Now, I want to minimize the RMSE between my model and the experimental values by changing m.d value which is the fitting parameter in the pde.

I incorporated the relevant code below.

--------

m.d = Var(initialize =50, bounds = (0.01, 100))

m.p = Var(m.a, m.x)

m.pi = Param(initialize=3.1416)

m.a = ContinuousSet(bounds=(-m.pi/2, m.pi/2))

m.x = ContinuousSet(bounds = (0, 1))

m.dpdx = DerivativeVar(m.p, wrt = m.x)

m.dpda2 = DerivativeVar(m.p, wrt = (m.a, m.a))

m.dpda = DerivativeVar(m.p, wrt = m.a)

def _pde(m, i, j):

    if j ==0: 

        return Constraint.Skip

    return m.ux[j]*m.dpdx[i,j] == m.d*m.dpda2[i,j] + m.r*(2*m.p[i,j]*cos(2*m.aa[i])+sin(2*m.aa[i])*m.dpda[i,j])

discretizer = TransformationFactory('dae.finite_difference')

discretizer2 = TransformationFactory('dae.collocation')

discretizer.apply_to(m, nfe=25, wrt=m.x, scheme='BACKWARD')

discretizer2.apply_to(m, nfe=170, ncp=3, wrt=m.a) 

def _obj(m):

    AAA = 0

    for i in range(17):    

        AAA = AAA + (sqrt(mean_squared_error([value(m.p[sorted(m.a)[10*i], sorted(m.x)[8]])], [YY1[i]])) +

                 sqrt(mean_squared_error([value(m.p[sorted(m.a)[10*i], sorted(m.x)[12]])], [YY2[i]])) +

                 sqrt(mean_squared_error([value(m.p[sorted(m.a)[10*i], sorted(m.x)[16]])], [YY3[i]])) +

                 sqrt(mean_squared_error([value(m.p[sorted(m.a)[10*i], sorted(m.x)[20]])], [YY4[i]])))

    return AAA

m.obj = Objective(rule = _obj, sense=minimize)

--------

Here, YY is the experimental value that I inserted as a numpy array.

My question is that:

If I want to find the best Dp value (m.d) that minimizes AAA, do I have to incorporate m.d like m.p = Var(m.a, m.x, m.d) in the above code?

Or by declaring m.d alone (m.d = Var(initialize =50, bounds = (0.01, 100))) is enough to do the optimization to find best m.d?

Thank you.

Jinhyeun Kim

[Nicholson, Bethany L.]

Do you intend for `m.d` to be a single value or to have different values depending on the index?

 

Bethany

--
You received this message because you are subscribed to the Google Groups "Pyomo Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pyomo-forum/9641a9c0-d326-437c-8772-2b57ded32505%40googlegroups.com.

[Jinhyeun Kim]

I intend m.d to be a single value for all index, but make it as a fitting parameter that I can optimize RMSE by changing m.d.

(In the PDE, I have been changing the m.d value and achieves optimum value that minimize RMSE between the experimental values.)

Below is the code that I am using now, and my code sometimes converge to local infeasibility, or it reached maximum iteration. (and the plot does not match with experimental data).

I was still figuring out the reason for infeasibility, and I was guessing the several possible reasons:

1. Discretization Method Problem. PDE I am using is the convection-dominated problem. Streamline Upwind Petrov-Galerkin (SUPG) finite element method is a popular method for solving this kind of PDE that eliminates the wiggling problem, but I want to know whether the collocation method also can solve this convection-dominated PDE.

2. Boundary Conditions Problem. I have been using six different boundary conditions (con1, con2, con3, con4, con5, _int) and combine them differently (con1+con2 / con1+con3 / con1 +con2+con4+con5+_int) for every run to find out the reason. My purpose is to reproduce the PDE in the paper, and I tested the boundary condition (con1 + con2) they used to discretize the equation (they used SUPG for discretization), which results in infeasibility for my code. So I am trying to test a different combination of the boundary condition which makes sense physically.

3. Solver Problem. I have been using ipopt solver with different linear_solver (what I used now is pardiso), and used BARON with the demo version which results in a license issue. 

Or I think it is just a simple coding mistake, maybe in my code I use different Dp values depending on the Index value. If I declared the m.d like this (m.d = Var(initialize = 2, bounds = (0.01, 100))), does it say m.d would vary depending on the indexed value (for example, are m.d different when m.a = 0 and m.a = 0.5)? then, how can I declare the m.d value that is used for optimization but not vary depending on the index?

Thank you for your help.

Jinhyeun Kim

 

m = ConcreteModel()

m.r = Param(initialize = 10) 

m.d = Var(initialize = 2, bounds = (0.01, 100))

m.pi = Param(initialize=3.1416)

m.a = ContinuousSet(bounds=(-m.pi/2, m.pi/2)) 

m.x = ContinuousSet(bounds = (0, 1)) 

m.p = Var(m.a, m.x)

m.dpdx = DerivativeVar(m.p, wrt = m.x)

m.dpda2 = DerivativeVar(m.p, wrt = (m.a, m.a))

m.dpda = DerivativeVar(m.p, wrt = m.a)

def _ux_rule(m, j):

    return 1/(1-(1-1/m.r)*(j))

m.ux = Expression(m.x, rule=_ux_rule)

def _a_rule(m, i):

    return i

m.aa = Expression(m.a, rule=_a_rule)

def _dudx_rule(m, j): 

    return (1-1/m.r)/((1-(1-1/m.r)*(j))**2)

m.dudx = Expression(m.x, rule=_dudx_rule)

def _pde(m, i, j):

    if j ==0: #or i == m.pi/2 or i == -m.pi/2:

        return Constraint.Skip

    return m.ux[j]*m.dpdx[i,j] == m.d*m.dpda2[i,j] + m.r*(2*m.p[i,j]*cos(2*m.aa[i])+sin(2*m.aa[i])*m.dpda[i,j])

m.pde = Constraint(m.a, m.x, rule=_pde)

def _con1(m, j):

    if j == 0:

        return Constraint.Skip

    return m.dpda[-m.pi/2, j] == m.dpda[m.pi/2, j]

m.con1 = Constraint(m.x, rule=_con1)

def _con2(m, i):

    return m.p[i, 0] == 1/m.pi

m.con2 = Constraint(m.a, rule=_con2)

def _con3(m, j):

    if j == 0:

        return Constraint.Skip

    return m.p[-m.pi/2, j] == m.p[m.pi/2, j]

m.con3 = Constraint(m.x, rule=_con3)

discretizer = TransformationFactory('dae.finite_difference')

discretizer2 = TransformationFactory('dae.collocation')

discretizer.apply_to(m, nfe=20, wrt=m.x, scheme='BACKWARD')

discretizer2.apply_to(m, nfe=100, ncp=10, wrt=m.a) # change nfe to 170

'''def _intp(m, i, j):

    return m.p[i, j]

m.intp = IndexedIntegral(m.a, m.x, wrt=m.a, rule=_intp)

def _con4(m, j):

    return m.intp[j] == 1

m.con4 = Constraint(m.x, rule=_con4)

def _con5(m, i, j):

    return m.p[i, j] >= 0

m.con5 = Constraint(m.a, m.x, rule=_con5)'''

def _obj(m): 

    cxz = {}

    for i in range(len(sorted(m.a))):

        cxz[(sorted(m.a)[i])] = est.predict(np.array(m.a).reshape(-1,1))[i][0]

        

    return sum((m.p[i, 0.75] - cxz[i])**2 for i in m.a) 

m.obj = Objective(rule = _obj, sense=minimize)

solver = SolverFactory('ipopt')

results = solver.solve(m, tee=True)

'''solver = SolverFactory('gurobi', solver_io='direct')

solver.solve(m)'''

#solver=SolverFactory('gams')

#solver.solve(m,solver='baron', tee=True)

#solver.solve(m,solver='XPRESS', tee=True)

print(value(m.d))