The story behind the tutorial model is as follows: A firm produces a single output using a production technology that combines labor, capital, and an aggregate material input according to a CES (constant elasticity of substitution) production function, i.e.

\[ Out = (A_{l} L^{-\rho} + A_{k} K^{-\rho} + A_{Inp} Inp^{-\rho})^{-1/\rho} \]

where:

\(Out\) is the volume of output,
\(L\) is the volume of labor input,
\(K\) is the volume of capital input,
\(Inp\) is the volume of material input, and
\(A_{l}, A_{k}, A_{Inp}\), and \(Rho\) are known parameters.

The demand for the output depends on the price according to the demand function

\[ Out = ( P_{0} - P ) / D \]

where:

\(P\) is the price, and
\(P_{0}\) and \(D\) are known parameters.

The firm is a family firm, and it has a certain amount of labor, \(L_{0}\), available from the family at no direct charge. Additional labor can be hired in any amount at a fixed wage rate, \(W\). The firm owns a certain amount of capital, \(K_{0}\), and it cannot change this in the short run. Material inputs can be purchased at a fixed unit price, \(P_{Inp}\), independent of quantity.

The objective of the firm is to maximize the short run profit, i.e. maximize sales revenues minus variable production costs, subject to the constraints and bounds described above. In the following we will use the term "bound" for simple lower and upper bounds on a single variable, and the term "constraint" for all other types of relationships involving nonlinearities and/or several variables, including the objective function. Although simple lower and upper bounds on single variables can be formulated as general inequalities, it is much more efficient to formulate them as “bounds” in a solver like CONOPT.

Mathematically, the optimization problem can be formulated as:

\[ \begin{array}[t]{rl} \max & \displaystyle Out P - (L-L_{0})W - Inp P_{Inp} \\ & \\ \text{subject to} & \displaystyle Out = (A_{l} L^{-\rho} + A_{k} K^{-\rho} + A_{Inp} Inp^{-\rho})^{-1/\rho} \\ & \displaystyle Out = ( P_{0} - P ) / D \\ & \displaystyle L \geq L_{0},\ K = K_{0},\ P \geq 0,\ Out \geq 0, \text{ and } Inp \geq 0 \end{array} \]

For the numerical example, we assume the following parameter values: \(W = 1.0, L_{0} = 0.1, P_{Inp} = 1.0, A_{l} = 0.16, A_{k} = 2.00, A_{Inp} = 0.16, Rho = 1.0, K_{0} = 4.0, P_{0} = 2.0,\) and \(D = 0.5\). \(Rho = 1.0\) corresponds to an elasticity of substitution between inputs of 0.5.

Using the above parameter values, the optimization problem is given by:

\[ \begin{array}[t]{rl} \max & \displaystyle Out P - (L-0.1) - Inp \\ & \\ \text{subject to} & \displaystyle Out = (0.16 L^{-1} + 0.5 + 0.16 Inp^{-1})^{-1} \\ & \displaystyle Out = ( 2 - P ) / 0.5 \\ & \displaystyle L \geq 0.1,\ P \geq 0,\ Out \geq 0, \text{ and } Inp \geq 0 \end{array} \]

Constraint reorganizations

There are four main steps that need to be followed when implementing a model for CONOPT. The first three are related to the specification and reorganization of constraints. The final step is related to the specification of the Jacobian, which will be discussed in the next section.

Reformulate the constraints

As a first step, the constraints need to be reformulated into the CONOPT format. First, move the variable terms to the left-hand side and the constant terms to the right hand side. In other words, rewrite each constraint in the form

\[ g(x) =\!\!R\!\!= b \]

where:

\(x\) represents the vector of optimization variables, in the following called model variables,
\(g(x)\) is a function of x,
\(b\) is a constant, and
\(=\!\!R\!\!=\) represent one of the four relationships: equal ( \(=\)), less than or equal ( \(\leq\)), greater than or equal ( \(\geq\)), or non-constraining ( \(=\!\!N\!\!=\)). Non-constraining is used for the objective function. It can also be used for alternative objective functions, or it can be used to turn constraints off without actually removing them.

Applied to the base model we get

\begin{align} &Out P - L W - Inp P_{Inp} &&=\!\!N\!\!= &&-L_{0}W \label{eqn:objective} \\ &(A_{l} L^{-\rho} + A_{k} K^{-\rho} + A_{Inp} Inp^{-\rho})^{-1/\rho} - Out &&= &&0 \label{eqn:constraint1} \\ &Out + P/D &&= &&P_{0}/D \label{eqn:constraint2} \end{align}

We have numbered the constraints at the left and we will refer to these numbers in the following. The constant term in the objective function, \(-L_{0}W\), can be ignored in an optimization. However, we keep the term in the model because it will allow us to read the real value of the objective function directly from the output produced by CONOPT. A constant term in a nonlinear constraint can also be included as part of the nonlinear function g(x); the right hand side will in this case become zero.

Identify linear and nonlinear terms

The second step in implementing an optimization model in CONOPT is to separate each individual left hand side into a sum of linear term and a remaining nonlinear term. Applied to our model we get:

For the objective function \(\eqref{eqn:objective}\)
- linear terms: \(- L W - Inp P_{Inp}\),
- nonlinear terms: \(Out P\).
For constraint \(\eqref{eqn:constraint1}\)
- linear terms: \(-Out\),
- nonlinear terms: \((A_{l} L^{-\rho} + A_{k} K^{-\rho} + A_{Inp} Inp^{-\rho})^{-1/\rho}\).
For constraint \(\eqref{eqn:constraint2}\)
- linear terms: \(Out + P/D\),
- nonlinear terms: none.
NOTE:** CONOPT uses the convention, that if a variable appears both linearly and non-linearly in a constraint expression, then the corresponding linear term must be included in the overall nonlinear term. Consider for example the expression \(e^{x} + 3x + 2y\). \(x\) appears nonlinearly and the linear term \(3x\) is therefore for the purpose of CONOPT counted as part of the overall nonlinear term together with \(e^{x}\); the only linear term is \(2y\). CONOPT also uses the convention that the separation into linear and nonlinear terms is done on an equation by equation basis. Therefore, the content of one equation does not influence the separation of the other equations.

Reference the constraints

The third step is to assign sequence numbers starting at 1 to the variables and constraints so they can be referenced within a Fortran framework using one-dimensional array notation. The constraints have been numbered above. For the variables we will use X[1] = \(L\), X[2] = \(Inp\), X[3] = \(Out\), and X[4] = \(P\). We have not included \(K\) as a model variable. The value of \(K\) is fixed, and in the following, we will simply consider \(K\) to be a parameter. It will make the model slightly faster to solve. On the other hand, we will not get any marginal cost (shadow price) information about \(K\) in the final solution.

For the constraints, these must also be numbered from 1. As such, the objective \(\eqref{eqn:objective}\) is indexed as 1, constraint \(\eqref{eqn:constraint1}\) is 2 and \(\eqref{eqn:constraint2}\) is 3.

The Jacobian Structure

The fourth step in implementing the model is to create a matrix that shows the structure of the matrix of derivatives, also known as the Jacobian of the model. The matrix must show which model variables appear in which constraints, and which relationships are linear and which are nonlinear. For the linear relationships, it must show the coefficient. In the below example the nonlinear cells are marked \(NL\). Many cells in the Jacobian will be empty. The matrix for our example is given by

\begin{equation*} \left[ \begin{array}{cccc} -1 & -1 & NL & NL \\ NL & NL & -1 & \\ & & +1 & +2 \\ \end{array} \right] \end{equation*}

There is one column for each model variable, ordered as \((L, Inp, Out, P)\) and one row for each constraint. Note, that the linearity/nonlinearity of the elements of the Jacobian reflects the separation of the constraint functions into linear and nonlinear terms as described above. The non-empty cells must be numbered column-wise starting from 0 for use in a sparse matrix storage. The vector of the non-empty cells is given by

\begin{equation*} [ \begin{array}{ccccccccc} -1 & NL & -1 & NL & NL & -1 & +1 & NL & +2 \end{array} ] \end{equation*}