The Xpress NonLinear API Functions
Instead of writing an extended MPS file and reading in the model from the file, it is possible to embed Xpress NonLinear directly into your application, and to create the problem, solve it and analyze the solution entirely by using the Xpress NonLinear API functions. This example uses the C header files and API calls. We shall assume you have some familiarity with the Xpress Optimizer API functions in XPRS.DLL.
The structure of the model and the naming system will follow that used in the previous section, so you should read the chapter Modeling in Extended MPS Format first.
Header files
The header file containing the Xpress NonLinear definitions is xslp.h. This must be included together with the Xpress Optimizer header xprs.h. xprs.h must come first.
#include "xprs.h" #include "xslp.h"
Initialization
Xpress NonLinear and Xpress Optimizer both need to be initialized, and an empty problem created. All Xpress NonLinear functions return a code indicating whether the function completed successfully. A non-zero value indicates an error. For ease of reading, we have for the most part omitted the tests on the return codes, but a well-written program should always test the values.
XPRSprob mprob; XSLPprob sprob; if (ReturnValue=XPRSinit(NULL)) goto ErrorReturn; if (ReturnValue=XSLPinit()) goto ErrorReturn; if (ReturnValue=XPRScreateprob(&mprob)) goto ErrorReturn; if (ReturnValue=XSLPcreateprob(&sprob, &mprob)) goto ErrorReturn;
Callbacks
It is good practice to set up at least a message callback, so that any messages produced by the system appear on the screen or in a file. The XSLPsetcbmessage function sets both the Xpress NonLinear and Xpress Optimizer callbacks, so that all messages appear in the same place.
XSLPsetcbmessage(sprob, XSLPMessage, NULL); |
void XPRS_CC XSLPMessage(XSLPprob my_prob, void *my_object, char *msg, int len, int msg_type) { switch (msg_type) { case 4: /* error */ case 3: /* warning */ case 2: /* dialogue */ case 1: /* information */ printf("%s\n", msg); break; default: /* exiting */ fflush(stdout); break; } }
This is a simple callback routine, which prints any message to standard output.
Creating the linear part of the problem
The linear part of the problem, and the definitions of the rows and columns of the problem are carried out using the normal Xpress Optimizer functions.
#define MAXROW 20 #define MAXCOL 20 #define MAXELT 50 int nRow, nCol, nSide, nRowName, nColName; int Sin, Cos; char RowType[MAXROW]; double RHS[MAXROW], OBJ[MAXCOL], Element[MAXELT]; double Lower[MAXCOL], Upper[MAXCOL]; int ColStart[MAXCOL+1], RowIndex[MAXELT]; char RowNames[500], ColNames[500];
In this example, we have set the dimensions by using #define statements, rather than working out the actual sizes required from the number of sides and then allocating the space dynamically.
nSide = 5; nRowName = 0; nColName = 0;
By making the number of sides a variable (nSide) we can create other polygons by changing its value.
It is useful – at least while building a model – to be able to see what has been created. We will therefore create meaningful names for the rows and columns. nRowName and nColName count along the character buffers RowNames and ColNames.
nRow = nSide-2 + (nSide-1)*(nSide-2)/2 + 1; nCol = (nSide-1)*2 + 2; for (i=0; i<nRow; i++) RHS[i] = 0;
The number of constraints is:
nSide-2 | for the relationships between adjacent thetas. |
(nSide-1)*(nSide-2)/2 | for the distances between pairs of vertices. |
1 | for the OBJEQ non-linear "objective function". |
The number of columns is:
nSide-1 | for the thetas. |
nSide-1 | for the rhos. |
1 | for the OBJX objective function column. |
1 | for the "equals column". |
We are using "C"-style numbering for rows and columns, so the counting starts from zero.
nRow = 0; RowType[nRow++] = 'E'; /* OBJEQ */ nRowName = nRowName + 1 + sprintf(&RowNames[nRowName], "OBJEQ"); for (i=1; i<nSide-1; i++) { RowType[nRow++] = 'G'; /* T2T1 .. T4T3 */ RHS[i] = 0.001; nRowName = nRowName + 1 + sprintf(&RowNames[nRowName], "T%dT%d", i+1, i); }
This sets the row type indicator for OBJEQ and the theta relationships, with a right hand side of 0.001. We also create row names in the RowNames buffer. Each name is terminated by a NULL character (automatically placed there by the sprintf function). sprintf returns the length of the string written, excluding the terminating NULL character.
for (i=1; i<nSide-1; i++) { for (j=i+1; j<nSide; j++) { RowType[nRow] = 'L'; RHS[nRow++] = 1.0; nRowName = nRowName + 1 + sprintf(&RowNames[nRowName], "V%dV%d", i, j); } }
This defines the L-type rows which constrain the distances between pairs of vertices. The right hand side is 1.0 (the maximum value) and the names are of the form ViVj.
for (i=0; i<nCol; i++) { OBJ[i] = 0; /* objective function */ Lower[i] = 0; /* lower bound normally zero */ Upper[i] = XPRS_PLUSINFINITY; /* upper bound = infinity */ }
This sets up the standard column data, with objective function entries of zero, and default bounds of zero to plus infinity. We shall change these for the individual items as required.
nCol = 0; nElement = 0; ColStart[nCol] = nElement; OBJ[nCol] = 1.0; Lower[nCol++] = XPRS_MINUSINFINITY; /* free column */ Element[nElement] = -1.0; RowIndex[nElement++] = 0; nColName = nColName + 1 + sprintf(&ColNames[nColName], "OBJX");
This starts the construction of the matrix elements. nElement counts through the Element and RowIndex arrays, nCol counts through the ColStart, OBJ, Lower and Upper arrays. The first column, OBJX, has the objective function value of +1 and a value of -1 in the OBJEQ row. It is also defined to be "free", by making its lower bound equal to minus infinity.
iRow = 0 for (i=1; i<nSide; i++) { nColName = nColName + 1 + sprintf(&ColNames[nColName], "THETA%d", i); ColStart[nCol++] = nElement; if (i < nSide-1) { Element[nElement] = -1; RowIndex[nElement++] = iRow+1; } if (i > 1) { Element[nElement] = 1; RowIndex[nElement++] = iRow; } iRow++; }
This creates the relationships between adjacent thetas. The tests on i are to deal with the first and last thetas which do not have relationships with both their predecessor and successor.
Upper[nCol-1] = 3.1415926;
This sets the bound on the final theta to be π. The column index is nCol-1 because nCol has already been incremented.
nColName = nColName + 1 + sprintf(&ColNames[nColName], "="); ColStart[nCol] = nElement; Lower[nCol] = Upper[nCol] = 1.0; /* fixed at 1.0 */ nCol++;
This creates the "equals column" – its name is "=" and it is fixed at a value of 1.0.
for (i=1; i<nSide; i++) { Lower[nCol] = 0.01; /* lower bound */ Upper[nCol] = 1; ColStart[nCol++] = nElement; nColName = nColName + 1 + sprintf(&ColNames[nColName], "RHO%d", i); } ColStart[nCol] = nElement;
The remaining columns – the rho variables – have only non-linear coefficients and so they do not appear in the linear section except as empty columns. They are bounded between 0.01 and 1.0 but have no entries. The final entry in ColStart is one after the end of the last column.
XPRSsetintcontrol(mprob, XPRS_MPSNAMELENGTH, 16);
If you are creating your own names – as we are here – then you need to make sure that Xpress Optimizer can handle both the names you have created and the names that will be created by Xpress NonLinear. Typically, Xpress NonLinear will create names which are three characters longer than the names you have used. If the longest name would be more than 8 characters, you should set the Xpress Optimizer name length to be larger – it comes in multiples of 8, so we have used 16 here. If you do not make the name length sufficiently large, then the XPRSaddnames function will return an error either here or during the Xpress NonLinear "construct" phase.
XPRSloadlp(mprob, "Polygon", nCol, nRow, RowType, RHS, NULL, OBJ, ColStart, NULL, RowIndex, Element, Lower, Upper);
This actually loads the model into Xpress Optimizer. We are not using ranges or column element counts, which is why the two arguments are NULL.
XPRSaddnames(mprob, 1, RowNames, 0, nRow-1); XPRSaddnames(mprob, 2, ColNames, 0, nCol-1);
The row and column names can now be added.
Adding the non-linear part of the problem
Be warned – this section is complicated, but it is the most efficient way – from SLP's point of view – to input formulae. See the next section for a much easier (but less efficient) way of inputting the formulae directly.
#define MAXTOKEN 200 #define MAXCOEF 20 ... int Sin, Cos; ColIndex[MAXCOL]; FormulaStart[MAXCOEF]; Type[MAXTOKEN]; double Value[MAXTOKEN], Factor[MAXCOEF];
The arrays for the non-linear part can often be re-used from the linear part. The new arrays are ColIndex (for the column index of the coefficients), FormulaStart and Factor for the coefficients, and Type and Value to hold the internal forms of the formulae.
XSLPgetindex(sprob, XSLP_INTERNALFUNCNAMES, "SIN", &Sin); XSLPgetindex(sprob, XSLP_INTERNALFUNCNAMES, "COS", &Cos);
We will be using the Xpress NonLinear internal functions SIN and COS. The XSLPgetindex function finds the index of an Xpress NonLinear entity (character variable, internal or user function).
nToken = 0; nCoef = 0; RowIndex[nCoef] = 0; ColIndex[nCoef] = nSide; Factor[nCoef] = 0.5; FormulaStart[nCoef++] = nToken;
For each coefficient, the following information is required:
RowIndex | the index of the row. |
ColIndex | the index of the column. |
FormulaStart | the beginning of the internal formula array for the coefficient. |
Factor | this is optional. If used, it holds a constant multiplier for the formula. This is particularly useful where the same formula appears in several coefficients, but with different signs or scaling. The formula can be used once, with different factors. |
for (i=1; i<nSide-1; i++) { Type[nToken] = XSLP_COL; Value[nToken++] = nSide+i+1; Type[nToken] = XSLP_COL; Value[nToken++] = nSide+i; Type[nToken] = XSLP_OP; Value[nToken++] = XSLP_MULTIPLY; Type[nToken] = XSLP_RB; Value[nToken++] = 0; Type[nToken] = XSLP_COL; Value[nToken++] = i+1; Type[nToken] = XSLP_COL; Value[nToken++] = i; Type[nToken] = XSLP_OP; Value[nToken++] = XSLP_MINUS; Type[nToken] = XSLP_IFUN; Value[nToken++] = Sin; Type[nToken] = XSLP_OP Value[nToken++] = XSLP_MULTIPLY; if (i>1) { Type[nToken] = XSLP_OP; Value[nToken++] = XSLP_PLUS; } }
This looks very complicated, but it is really just rather large. We are using the "reverse Polish" or "parsed" form of the formula for area. The original formula, written in the normal way, would look like this:
RHO2 * RHO1 * SIN ( THETA2 - THETA1 ) + .......
In reverse Polish notation, tokens are pushed onto the stack or popped from it. Typically, this means that a binary operation A x B is written as A B x (push A, push B, pop A and B and push the result). The first term of our area formula then becomes:
RHO2 RHO1 * ) THETA2 THETA1 - SIN *
Notice that the right hand bracket appears as an explicit token. This allows the SIN function to identify where its argument list starts – and incidentally allows functions to have varying numbers of arguments.
Each token of the formula is written as two items – Type and Value.
Type is an integer and is one of the defined types of token, as given in the xslp.h header file. XSLP_CON, for example, is a constant; XSLP_COL is a column.
Value is a double precision value, and its meaning depends on the corresponding Type. For a Type of XSLP_CON, Value is the constant value; for XSLP_COL, Value is the column number; for XSLP_OP (arithmetic operation), Value is the operand number as defined in xslp.h; for a function (type XSLP_IFUN for internal functions, XSLP_FUN for user functions), Value is the function number.
A list of tokens for a formula is always terminated by a token of type XSLP_EOF.
The loop writes each term in order, and adds terms (using the XSLP_PLUS operator) after the first pass through the loop.
for (i=1; i<nSide-1; i++) { for (j=i+1; j<nSide; j++) { RowIndex[nCoef] = iRow++; ColIndex[nCoef] = nSide; Factor[nCoef] = 1.0; FormulaStart[nCoef++] = nToken; Type[nToken] = XSLP_COL; Value[nToken++] = nSide+i; Type[nToken] = XSLP_CON; Value[nToken++] = 2; Type[nToken] = XSLP_OP; Value[nToken++] = XSLP_EXPONENT; Type[nToken] = XSLP_COL; Value[nToken++] = nSide+j; Type[nToken] = XSLP_CON; Value[nToken++] = 2; Type[nToken] = XSLP_OP; Value[nToken++] = XSLP_PLUS; Type[nToken] = XSLP_CON; Value[nToken++] = 2; Type[nToken] = XSLP_COL; Value[nToken++] = nSide+i; Type[nToken] = XSLP_OP; Value[nToken++] = XSLP_MULTIPLY; Type[nToken] = XSLP_COL; Value[nToken++] = nSide+j; Type[nToken] = XSLP_OP; Value[nToken++] = XSLP_MULTIPLY; Type[nToken] = XSLP_RB; Value[nToken++] = 0; Type[nToken] = XSLP_COL; Value[nToken++] = j; Type[nToken] = XSLP_COL; Value[nToken++] = i; Type[nToken] = XSLP_OP; Value[nToken++] = XSLP_MINUS; Type[nToken] = XSLP_IFUN; Value[nToken++] = Cos; Type[nToken] = XSLP_OP Value[nToken++] = XSLP_MULTIPLY; Type[nToken] = XSLP_OP; Value[nToken++] = XSLP_MINUS; Type[nToken] = XSLP_EOF; Value[nToken++] = 0; } }
This writes the formula for the distances between pairs of vertices. It follows the same principle as the previous formula, writing the formula in parsed form as:
RHOi 2 RHOj 2 + 2 RHOi * RHOj * ) THETAj THETAi - COS * -
XSLPloadcoefs(sprob, nCoef, RowIndex, ColIndex, Factor, FormulaStart, 1, Type, Value); |
The XSLPloadcoefs is the most efficient way of loading non-linear coefficients into a problem. There is an XSLPaddcoefs function which is identical except that it does not delete any existing coefficients first. There is also an XSLPchgcoef function, which can be used to change individual coefficients one at a time. Because we are using internal parsed format, the "Parsed" flag in the argument list is set to 1.
Adding the non-linear part of the problem using character formulae
Provided that all entities – in particular columns and user functions – have explicit and unique names, the non-linear part can be input by writing the formulae as character strings. This is not as efficient as using the XSLPloadcoefs() function but is generally easier to understand.
/* Build up nonlinear coefficients */ /* Allow space for largest formula - approx 50 characters per side for area */ CoefBuffer = (char *) malloc(50*nSide);
We shall be using large formulae, so we need a character buffer large enough to hold the largest formula we are using. The estimate here is 50 characters per side of the polygon for the area formula, which is the largest we are using.
/* Area */ Factor = 0.5; BufferPos = 0; for (i=1; i<nSide-1; i++) { if (i > 1) { BufferPos = BufferPos + sprintf(&CoefBuffer[BufferPos], " + "); } BufferPos = BufferPos + sprintf(&CoefBuffer[BufferPos], "RHO%d * RHO%d * SIN ( THETA%d - THETA%d )", i+1, i, i+1, i); } XSLPchgccoef(sprob, 0, nSide, &Factor, CoefBuffer);
The area formula is of the form:
(RHO2*RHO1*SIN(THETA2-THETA1) + RHO3*RHO2*SIN(THETA3-THETA2) + ... ) / 2
The loop writes the product for each consecutive pair of vertices and also puts in the "+" sign after the first one.
The XSLPchgccoef function is a variation of XSLPchgcoef but uses a character string for the formula instead of passing it as arrays of tokens. The arguments to the function are:
RowIndex | the index of the row. |
ColIndex | the index of the column. |
Factor | this is optional. If used, it holds the address of a constant multiplier for the formula. This is particularly useful where the same formula appears in several coefficients, but with different signs or scaling. The formula can be used once, but with different factors. To omit it, use a NULL argument. |
CoefBuffer | the formula, written in character form. |
In this case, RowIndex is zero and ColIndex is nSide (the "equals" column).
/* Distances */ Factor = 1.0; for (i=1; i<nSide-1; i++) { for (j=i+1; j<nSide; j++) { sprintf(CoefBuffer, "RHO%d ^ 2 + RHO%d ^ 2 - 2 * RHO%d * RHO%d * COS ( THETA%d - THETA%d )", j, i, j, i, j, i); XSLPchgccoef(sprob, iRow, nSide, &Factor, CoefBuffer); iRow++; }
This creates the formula for the distance between pairs of vertices and writes each into a new row in the "equals" column.
Provided you have given names to any user functions in your program, you can use them in a formula in exactly the same way as SIN and COS have been used above.
Checking the data
Xpress NonLinear includes the function XSLPwriteprob which writes out a non-linear problem in text form which can then be checked manually. Indeed, the problem can then be run using the XSLP console program, provided there are no user functions which refer back into your compiled program. In particular, this facility does allow small versions of a problem to be checked before moving on to the full size ones.
XSLPwriteprob(sprob, "testmat", ""); |
The first argument is the Xpress NonLinear problem pointer; the second is the name of the matrix to be produced (the suffix ".mat" will be added automatically). The last argument allows various different types of output including "scrambled" names – that is, internally-generated names will be used rather than those you have provided. For checking purposes, this is obviously not a good idea.
Solving and printing the solution
XSLPmaxim(sprob, ""); |
The XSLPmaxim and XSLPminim functions perform a non-linear maximization or minimization on the current problem. The second argument can be used to pass flags as defined in the Xpress NonLinear Reference Manual.
XPRSwriteprtsol(mprob); |
The standard Xpress Optimizer solution print can be obtained by using the XPRSwriteprtsol function. The row and column activities and dual values can be obtained using the XPRSgetsol function.
In addition, you can use the XSLPgetvar function to obtain the values of SLP variables – that is, of variables which are in non-linear coefficients, or which have non-linear coefficients. If you are using cascading (see the Xpress NonLinear reference manual for more details) so that Xpress NonLinear recalculates the values of the dependent SLP variables at each SLP iteration, then the value from XSLPgetvar will be the recalculated value, whereas the value from XPRSgetsol will be the value from the LP solution (before recalculation).
Closing the program
The XSLPdestroyprob function frees any system resources allocated by Xpress NonLinear for the specific problem. The problem pointer is then no longer valid. XPRSdestroyprob performs a similar function for the underlying linear problem mprob. The XSLPfree function frees any system resources allocated by Xpress NonLinear. You must then call XPRSfree to perform a similar operation for the optimizer.
XSLPdestroyprob(sprob); XPRSdestroyprob(mprob); XSLPfree(); XPRSfree();
If these functions are not called, the program may appear to have worked and terminated correctly. However, in such a case there may be areas of memory which are not returned to the system when the program terminates and so repeated executions of the program will result in progressive loss of available memory to the system, which will manifest iself in poorer performance and could ultimately produce a system crash.
Adding initial values
So far, Xpress NonLinear has started by using values which it estimates for itself. Because most of the variables are bounded, these initial values are fairly reasonable, and the model will solve. However, in general, you will need to provide initial values for at least some of the variables. Initial values, and other information for SLP variables, are provided using the XSLPloadvars function.
int VarType[MAXCOL]; double InitialValue[MAXCOL];
To load initial values using XSLPloadvars, we need an array (InitialValue) to hold the initial values, and a VarType array which is a bitmap to describe what information is being set for each variable.
for(i=1; i<nSide; i++) { ... InitialValue[nCol] = 3.14159*((double)i) / ((double)nSide); VarType[nCol] = 4; ... } ... for(i=1; i<nSide; i++) { InitialValue[nCol] = 1; VarType[nCol] = 4; }
These sections extend the loops for the columns in the earlier example. We set initial values for the thetas so that the vertices are spaced at equal angles; the rhos are all started at 1. We do not need to set a value for the equals column, because it is fixed at one. However, it is good practice to do so. In each case we set VarType to 4 because (as described in the Xpress NonLinear Reference Manual) Bit 2 of the type indicates that the initial value is being set.
for(i=0; i<nCol; i++) ColIndex[i] = i XSLPloadvars(sprob, nCol-1, &ColIndex[1], &VarType[1], NULL, NULL, NULL, &InitialValue[1], NULL);
XSLPloadvars can take several other arguments apart from the initial value. It is a general principle in Xpress NonLinear that using NULL for an argument means that there is no information being provided, and the current or default value will not be changed.
Because we built up the initial values as we went, the VarType and InitialValue arrays include column 0, which is OBJX and is not an SLP variable. As all the rest are SLP variables, we can simply start these arrays at the second item, and reduce the variable count by 1.
User functions
The most complicated formula in this model is the area calculation. With only 5 sides, it is still possible to write it out explicitly, but it becomes large (and perhaps inefficient) if the number of sides increases. The alternative is to calculate the formula in a function and then use the function within the model.
A user function is essentially a function which is not built in to Xpress NonLinear. It can be written in a language such as C or Fortran, and compiled into a DLL; it can be written as a set of formulae in an Excel spreadsheet (with or without a macro as well); it can be written entirely within an Excel macro. This example shows the area function written as a compiled C function.
A user function in C
This function calculates the area from an array of values, ordered as (RHO1, THETA1, RHO2, THETA2, ...). The number of items in the Values array is given as the first item in nArg.
double XPRS_CC MyFunc(double *Values, int *nArg) { int i; double Area; Area = 0; for(i=3; i<nArg[0]; i=i+2) { Area = Area + Values[i-3]*Values[i-1]*sin(Values[i]-Values[i-2]); } return Area*0.5; }
This is the standard interface for a user function in Xpress NonLinear. The first argument is an array of double precision values holding the values of the arguments for the Xpress NonLinear function in order; the second argument is an array of integers, the first of which contains the size of the first array.
The function must be declared using XPRS_CC as shown, to ensure that the correct function linkage is created.
This function can be compiled into a DLL. To make use of it, we also need to be able to access the formula from outside, so you may need to add suitable externalization definitions. In Visual C++ under Microsoft Windows, you can use a Definition File, containing an EXPORTS section, such as:
EXPORTS MyFunc=_MyFunc@8
Extending the polygon model
We can now declare this function in the model and use it instead of the explicit area formula.
nToken = 0; XSLPsetstring(sprob, &i, "MyFunc"); Type[nToken] = XSLP_STRING; Value[nToken++] = (double) i; Type[nToken] = XSLP_UFEXETYPE; Value[nToken++] = (double) 0x01; Type[nToken] = XSLP_UFARGTYPE; Value[nToken++] = (double) 023; XSLPsetstring(sprob, &i, "MyDLL.DLL"); Type[nToken] = XSLP_STRING; Value[nToken++] = (double) i; Type[nToken] = XSLP_EOF; XSLPloaduserfuncs(sprob, 1, Type, Value); XSLPaddnames(sprob, XSLP_USERFUNCNAMES, "MyArea", 1, 1);
User functions are declared using XSLPloaduserfuncs. The definition of the function is stored in parsed arrays similar to the ones used for defining formulae. There are two special token types used here; see the Xpress NonLinear Reference Manual for full details about the corresponding values.
We must also define the name of the function. This is a character string, and it is the first item in the array of tokens. To pass a character string to Xpress NonLinear, use the XSLPsetstring function to store the string and return an index to the string. Then use the index with the XSLP_STRING token type.
Because this is a DLL function, we must also define the name of the DLL. This is the first string after the tokens defining the function and argument types. For other types of function (for example, Excel spreadsheets or macros), other string parameters may be needed as well.
The XSLPaddnames function creates a name for the function to be used inside Xpress NonLinear when the function is referenced. It is what you will see if you write the problem out using XSLPwriteprob. It can be the same name as the function name in the DLL, but it does not have to be. If you are not writing the problem out, then you do not need to set a name at all.
Type[nToken] = XSLP_RB; Value[nToken++] = 0; for (i=nSide-1; i>0; i--) { Type[nToken] = XSLP_COL; Value[nToken++] = i; Type[nToken] = XSLP_COL; Value[nToken++] = nSide+i; } Type[nToken] = XSLP_FUN; Value[nToken++] = 1; Type[nToken] = XSLP_EOF; Value[nToken++] = 0;
In reverse Polish, the arguments to the function must appear in reverse order, so the items start with THETA4 and work down to RHO1. The arguments are preceded by a right bracket token and followed by the user function token for function number 1.
Internal user functions
The example above used a function written in a DLL. If the function is compiled into something else – for example, the main executable program – or is not externalized, then you will need to define its address explicitly.
void *Func; Func = MyFunc; XSLPchguserfuncaddress(sprob, 1, &Func);
XSLPchguserfuncaddress takes as its arguments the number of the function, and a pointer to its address. As usual, if the pointer is NULL, the data is left unaltered. The main use of the routine is to define the address of a user function directly, without relying on Xpress NonLinear to find it.