Pyomo dae for feasibility problem does not have same number of equations as unknowns

[Rob]

Hi all,

I've posted this question on Stackoverflow but have copied the post here:

I'm running the heat transfer example problem for pyomo.dae found in the pyomo.dae paper:

Nicholson, B., Siirola, J.D., Watson, JP. et al. Math. Prog. Comp. (2018) 10: 187. https://doi.org/10.1007/s12532-017-0127-0

I have been able to run the code correctly after copying and pasting it from the paper, but I'm confused about the results I'm getting back from ipopt

from pyomo.environ import *
from run_util import choose_executable  # my executable path on my machine
from pyomo.dae import *

m = ConcreteModel()
m.pi = Param(initialize=3.1416)
m.t = ContinuousSet(bounds=(0, 2))
m.x = ContinuousSet(bounds=(0, 1))
m.u = Var(m.x, m.t)
# Declare derivatives in the model
m.dudx = DerivativeVar(m.u, wrt=m.x)
m.dudx2 = DerivativeVar(m.u, wrt=(m.x, m.x))
m.dudt = DerivativeVar(m.u, wrt=m.t)

# Declare PDE
def _pde(m, i, j):
    if i == 0 or i == 1 or j == 0:
        return Constraint.Skip
    return m.pi ** 2 * m.dudt[i, j] == m.dudx2[i, j]
m.pde = Constraint(m.x, m.t, rule=_pde)

def _initcon(m, i):
    if i == 0 or i == 1:
        return Constraint.Skip
    return m.u[i, 0] == sin(m.pi * i)
m.initcon = Constraint(m.x, rule=_initcon)

def _lowerbound(m, j):
    return m.u[0, j] == 0
m.lowerbound = Constraint(m.t, rule=_lowerbound)

def _upperbound(m, j):
    return m.pi * exp(-j) + m.dudx[1, j] == 0
m.upperbound = Constraint(m.t, rule=_upperbound)
m.obj = Objective(expr=1)

# Discretize using Finite Difference and Collocation
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 = 20,ncp = 3,wrt = m.t)

solver = SolverFactory('ipopt', executable=choose_executable())
results = solver.solve(m, tee=True)

Since this is a fully specified problem, I expect the number of equations to be equal to the number of unknowns. However, the output from ipopt tells me that there are 60 more unknowns than equations. This doesn't seem correct to me. Is pyomo accounting for the derivative boundary condition correctly? Did I install ipopt incorrectly?

Below is the output I get from ipopt.

Ipopt 3.12: 
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************
This is Ipopt version 3.12, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...:    21257
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        0
Total number of variables............................:     6195
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:     6135
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.0000000e+00 3.14e+00 0.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.0000000e+00 1.68e-09 9.84e-04  -1.0 9.84e+00  -4.0 1.00e+00 1.00e+00h  1
   2  1.0000000e+00 1.10e-10 2.96e-11  -5.7 8.89e-07  -4.5 1.00e+00 1.00e+00h  1
Number of Iterations....: 2
                                   (scaled)                 (unscaled)
Objective...............:   1.0000000000000000e+00    1.0000000000000000e+00
Dual infeasibility......:   2.9639427023830905e-11    2.9639427023830905e-11
Constraint violation....:   1.0986800358381288e-10    1.0986800358381288e-10
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   1.0986800358381288e-10    1.0986800358381288e-10
Number of objective function evaluations             = 3
Number of objective gradient evaluations             = 3
Number of equality constraint evaluations            = 3
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 3
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 2