Optimization solver interfaces

[stroitzsch]

Optimization problem setup through cvxpy is slow for large test cases in FLEDGE (which seems to be a known issue: https://github.com/cvxpy/cvxpy/issues/708, https://github.com/cvxpy/cvxpy/issues/1159, https://github.com/cvxpy/cvxpy/issues/706). I'd like to try adding direct optimization solver interfaces to see if we can get better performance. To keep things clean & make switching interfaces easier in future, we will need some internal abstractions to define optimization variables / constraints / objective.

• [x] Refactor optimization problem standard form object with methods for variable / constraint / objective definition.
• [x] Refactor optimization problem standard form object to allow different solver interfaces (while retaining cvxpy as fallback interface).
• [x] Add gurobipy interface.
• [x] Enable retrieving of dual variables (requires labeling of constraints & retrieving dual values from solvers).
• [x] Enable adding fully featured matrix equations. (= broadcasting of multiplication along one dimension)
• [x] Enable quadratic objective (generate Q matrix).
• [x] Fix: Invalid constraint dimension if no constant defined.
• [x] Add high-level define_objective().
• [x] Enable using parameters in place of numerical values (For updating optimization problem without redefining constraints).
• [x] Replace optimization problem object with standard form object,
• [x] Redefine model class methods: define_optimization_variables(), define_optimization_constraints(), define_optimization_objective(), get_optimization_results(), get_optimization_dlmps().

[stroitzsch]

This is being implemented on branch feature/solver_interfaces and in pull request #12.

[tomschelo]

@stroitzsch I just pushed my first attempt, however, I cannot get it to work properly. When trying to solve, it returns "unfeasible or unbounded". However, Gurobi's "Coefficient statistics" are equal to the ones produced by the equivalent problem solved with cvxpy.

One key difference I encountered is that gurobipy forces your to define "x_vector" as a 1-D array --> shape=(len(variables), ) instead of shape=(len(variables), 1). Otherwise you cannot formulate a matrix-based constraint. Maybe I did something wrong in reshaping x and b?

[stroitzsch]

@tomschelo I think the issue was in the variable definition.

Looks like Gurobi defaults the lower bound to zero:

So we've got to explicitly set it to negative infinity:

x_vector = (
    gurobipy_problem.addMVar(
        shape=(len(self.variables), ),
        lb=-np.inf,
        ub=np.inf,
        vtype=GRB.CONTINUOUS,
        name='x_vector'
    )
)

I'll fix and push this.

[stroitzsch]

Retrieving duals & broadcasting done!

[kaiATtum]

whats broadcasting？

[stroitzsch]

@kaiATtum Check standard_form_example: Basically repeated matrix operation along an index dimension. E.g. state equation repeated for all timesteps. Broadcasting implementation avoids for-loops like for timestep in timesteps.

[stroitzsch]

@tomschelo I've fixed the grid constraint equations: Basically there was as in error in the way that constants were repeated for all timesteps with np.repeat(). To avoid this in future, I've enabled broadcasting for constants. Looks a bit hacky, because you have to set an index key set for the constants, but works fine.

I've also slightly rewritten the equations. The new formulation uses diagonal matrix multiplications to replace some of of the numpy broadcasting operations. This way we retain the sparse matrices, which may be useful for performance with bigger problems.

[stroitzsch]

I've also added "parameters" in the latest commits: Parameter values can be set with define_parameter(name, value), then the name string can be used in place of numerical values in define_constraint()

[stroitzsch]

Redefine model class methods also done. - Some testing / bugfixing / cleanup remains to be done.