BCLS Example Problem

A constrained linear least-squares problem and its formulation for BCLS, using the C language interface, are presented. This example is provided in MATLAB as an "M-file," bcls_ex1.m, in Fortran as bcls_ex1.for and in C as bcls_ex1.c on the BCLS_DLL distribution disk.

Problem Statement

A set of data points {x[i], y[i]}, i=1,…,m is approximated by a polynomial

p(x) = Sum{i=0, …,n}{a[i]*x[i]^(i)}.

The independent variable values x[I] satisfy 0<=x[I]<=1. Due to requirements of the model the value must satisfy p(0)=0 and additionally p'(0)=p'(1).

This is a classic example for linear least squares solving. The added twist is the pair of constraints. Once we have chosen a degree n=3,4,5, or 6, the linear problem is solved using the BCLS software.

C Program for BCLS Example Problem

#include <math.h>
#include <stdio.h>
#include "bclsc.h"

#define LDA     6
#define LDC     2
#define NDATA   6
#define NMAX    6

/*
        Fit data points by a polynomial.
     #0 Corresponds to a polynomial, with no constraints.
     #1 Constrains the polynomial to have value zero at zero.
     #2 Constraints as in #1 but with the periodic derivatives.

        Thus p(0)=0; also p'(0)=p'(1).
     This leads to three separate least squares problems with
     0, 1 and 2 constraints.  The degree varies from 3 to 6.
     This is Example 1 of WT, Inc. for the solver bclsc().
*/


void main() {
long int      i, j, k, n, ndeg, inform;
long int      id;
double        a[(NMAX+1)*LDA], c[(NMAX + 1)*LDC], coeff[NMAX + 1],
              val, xlb[NMAX + 1], xub[NMAX + 1], ylb[2], yub[2];
static double f[NDATA] = {0.2, 1.2214, 1.4918, 1.8221, 2.2255, 2.7183};
static double x[NDATA] = {0, .2, .4, .6, .8, 1.0};

    /* Define data to fit: */
    /* The degree of the polynomial: */
    for (ndeg = 3; ndeg <= 6; ndeg++)  {

        n = ndeg + 1;   /* The number of unknowns. */
        /* Fill in the matrix for the polynomial fitting: */
        for (i = 0; i < NDATA; i++)  {
            a[i*n]     = 1e0;
            a[i*n + 1] = x[i];
            for (j = 2; j < n; j++)  {
                a[i*n + j] = x[i]*a[i*n + j-1];
                }
            }

        /* Set no bounds for the coefficient sizes. */
        for (j = 0; j < n; j++)  {
            id = 2;
            xlb[j] = -DMACH (&id);
            xub[j] = DMACH (&id);
            }

        /* Fit data first with no constraints. */
        k = 0;  /* This is fit #0. */

        /* Compute the least-squares solution with no 
           constraints. */
        bclsc (NDATA, n, a, f, xlb, xub, k,
               c, ylb, yub, coeff, &inform);

        /* Add the constraint that p(0) = 0.  This can also be 
           achieved by a lower and upper bound of value 0 for 
           COEFF(1). */
        k = 1;  /* This is fit #1. */
        c[0] = 1e0;
        for (j = 1; j < n; j++)  {
            c[j] = 0;
            }
        ylb[0] = 0e0;
        yub[0] = 0e0;

        /* Compute the least-squares solution with ONE 
           constraint. */
        bclsc (NDATA, n, a, f, xlb, xub, k,
               c, ylb, yub, coeff, &inform);

        /* Add the constraint that p'(0) = p'(1).  This can also 
           be achieved by a lower and upper bound of value 0 for 
           COEFF(1). */
        k = 2;  /* This is fit #3. */
        c[n]     = 0e0;
        c[n + 1] = 0e0;
        for (j = 2; j < n; j++)  {
            c[n + j] = j;
            }
        ylb[1] = 0e0;
        yub[1] = 0e0;

        /* Compute the least-squares solution with TWO 
           constraints. */
        bclsc (NDATA, n, a, f, xlb, xub, k,
               c, ylb, yub, coeff, &inform);

        printf ("ndeg = %ld     inform = %ld\n", ndeg, inform);

        /* Check that p(0)=0 and p'(0)=p'(1)., approximately. */
        val = fabs (coeff[0]); /* This is |p(0)|. */
        for (j = 2; j < n; j++)  {
            val = val + j * coeff[j];
            }
        id = 4;

        if (fabs (val) > sqrt (DMACH (&id)))  {
            printf ("The conditions p(0)=0 and p'(0)=p'(1) "
                    "are violated.\n");
            }

    }
    printf ("BclsC_ex1 finished successfully.\n");
    MessageBox ((HWND) NULL, "BclsC_ex1 finished successfully",   
                             "BclsC Example Program", 
                              MB_OK+MB_ICONSTOP+MB_TASKMODAL );


}   /* end bclsc_ex1 */