-------- Forwarded Message --------
From: [email protected]
To: [email protected]
Subject: primal and dual simplex solves multiobjective problems
Date: Mon, 07 Dec 2015 15:34:33 +0100
Hi,
The attached patch against version 4.57 implements a kind of "multiobjective"
optimization. Given the "main" objective c_1,...,c_n and "auxiliary"
objectives <e_11,...e_1n>, <e_21,...,e_2n>,...<e_p1,...,e_pn>, find
the optimal value of the "main" objective, then within the solution
space minimize (maximize) the first auxiliary objective, within the
solution space minimize (maximize) the second auxiliary objective,
etc. It is a kind of generalization of asking the lexicographically
minimal optimal solution of the LP problem.
Please consider this patch as a tentative contribution to the package;
use it as you like. I would be happy to receive comments what I did
wrongly.
The following new api functions are added:
glp_set_multiobj_number(P,nobjs)
sets the number of auxiliary objectives; set to zero to delete
previous objectives.
glp_set_multiobj_coef(P,objno,j,coef)
set the objno-th auxiliary objective coefficient
glp_get_multiobj_coef(P,objno,j)
retrieve the coefficient
A new parameter in glp_smcp is parm->mobj, which should be set to
GLP_ON to do multiobjective optimization.
All added lines are between #ifdef CSL and #endif; there is only one
disabled line which checks the consistency of bound flags in the dual
problem.
In the primal/dual algorithm, primal_simplex()/dual_simplex() is
called repeatedly after calling next_objective() which goes over all
non-basic variables and fixes the value of those which have active
bound (by making the lower and upper bounds the same), and then
replaces the objective values with the next objective.
It is not clear whether csa->phase has to be set to zero in the dual
problem; and whether the when PSE pricing is on, whether it should be
invalidated from one invocation to the next.
Laszlo
Laszlo Csirmaz
CEU,
Alfred Renyi Mathematical Institute
Budapest
----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.
diff -r -c a/glpapi01.c b/glpapi01.c
*** a/glpapi01.c 2015-11-08 08:00:00.000000000 +0100
--- b/glpapi01.c 2015-12-03 14:31:51.078419642 +0100
***************
*** 83,88 ****
--- 83,95 ----
lp->row = xcalloc(1+lp->m_max, sizeof(GLPROW *));
lp->col = xcalloc(1+lp->n_max, sizeof(GLPCOL *));
lp->r_tree = lp->c_tree = NULL;
+ #ifdef CSL
+ /* multi objective part */
+ lp->cobj_num = 0;
+ lp->cobj_max = lp->cobj_end = 0;
+ lp->cobj_idx = NULL;
+ lp->cobj_val = NULL;
+ #endif
/* basis factorization */
lp->valid = 0;
lp->head = xcalloc(1+lp->m_max, sizeof(int));
***************
*** 221,226 ****
--- 228,275 ----
return;
}
+ #ifdef CSL
+ /***********************************************************************
+ * NAME
+ *
+ * glp_set_multiobj_number - set (change) number of additional objectives
+ *
+ * SYNOPSIS
+ *
+ * void glp_set_multiobj_number(glp_prob *lp, int nobjs);
+ *
+ * DESCRIPTION
+ *
+ * The routine glp_set_obj_number sets (changes) the number of additional
+ * objectives. When set to zero, no extra objectives are used, or the
+ * previously set objectives are deleted. New objectives are initially have
+ * all zero coefficients. */
+
+ void glp_set_multiobj_number(glp_prob *lp, int nobjs)
+ { glp_tree *tree = lp->tree;
+ int idx;
+ if (tree != NULL && tree->reason != 0)
+ xerror("glp_set_multiobj_number: operation not allowed\n");
+ if (nobjs < 0 )
+ xerror("glp_set_multiobj_number: nobjs=%d; invalid number\n", nobjs);
+ if (nobjs > M_MAX)
+ xerror("glp_set_multiobj_number: nobjs=%d; too many
objectives\n",nobjs);
+ if (nobjs == 0) /* delete all */
+ { lp->cobj_num = 0;
+ lp->cobj_end = 0;
+ }
+ else if (nobjs < lp->cobj_num) /* shrink */
+ { lp->cobj_num = nobjs;
+ for (idx=lp->cobj_idx[nobjs-1];lp->cobj_idx[idx]>=0;idx++);
+ lp->cobj_end = idx+1;
+ }
+ else if (nobjs > lp->cobj_num) /* expand */
+ { for (idx=lp->cobj_num+1;idx<=nobjs;idx++)
+ glp_set_multiobj_coef(lp,idx,0,0.0);
+ }
+ return;
+ }
+ #endif
/***********************************************************************
* NAME
*
***************
*** 691,696 ****
--- 740,855 ----
return;
}
+ #ifdef CSL
+ /***********************************************************************
+ * NAME
+ *
+ * glp_set_multiobj_coef - set (change) obj. coefficient or constant term
+ *
+ * SYNOPSIS
+ *
+ * void glp_set_multiobj_coef(glp_prob *lp, int objn, int j, double coef);
+ *
+ * DESCRIPTION
+ *
+ * The routine glp_set_multiobj_coef sets (changes) the objn'th objective
+ * coefficient at j-th column (structural variable) of the specified problem
+ * object.
+ *
+ * If the parameter j is 0, the routine sets (changes) the constant term
+ * ("shift") of the objective function. */
+
+
+ /* Non-zero coefficients of additional objectives are stored in
+ * cobj_idx[0:cobj_max-1] and cobj_val[0:cobj_max-1]. Objectives are indexed
+ * from 1 to obj_num. Coefficients for objective objno is stored as follows.
+ * cobj_val[objno-1] contains the "shift" (not used)
+ * cobj_idx[objno-1] contains a forward index "idx"
+ * starting from "idx" cobj_idx[idx] is the column number
+ * cobj_val[idx] is the objective value
+ * until cobj_idx[idx]==-1
+ *
+ * Alternative: cobj_val[objno-1] contains the objective value; the
+ * "shift" is stored at cobj_val[idx] where cobj_idx[idx]==-1.
+ *
+ * The routine mk_multiobjs_space(lp,idx) adds an empty slot at index idx
+ * by shifting cobj_idx and cobj_val by one to the right, and adjusting
+ * the pointers at 0:cobj_num-1. */
+ static void mk_multiobj_space(glp_prob *lp, int idx)
+ { int cend = lp->cobj_end; /* index of first free slot */
+ int i;
+ xassert(idx <= lp->cobj_end);
+ xassert(lp->cobj_num <= idx);
+ ++ lp-> cobj_end; /* now we have one more used elements */
+ if (cend >= lp->cobj_max)
+ /* no more space, expand both cobj_idx and cobj_val */
+ { lp->cobj_max += 200;
+ if (lp->cobj_max == 200)
+ { lp->cobj_idx = talloc(lp->cobj_max, int);
+ lp->cobj_val = talloc(lp->cobj_max, double);
+ }
+ else
+ { lp->cobj_idx = trealloc(lp->cobj_idx, lp->cobj_max, int);
+ lp->cobj_val = trealloc(lp->cobj_val, lp->cobj_max, double);
+ }
+ }
+ /* shift idx .. cend to the right by one */
+ for (i=cend;i>idx;i--)
+ { lp->cobj_idx[i]=lp->cobj_idx[i-1];
+ lp->cobj_val[i]=lp->cobj_val[i-1];
+ }
+ /* adjust pointers at the beginning */
+ for (i=0;i<lp->cobj_num;i++)
+ { if (lp->cobj_idx[i] >= idx)
+ lp->cobj_idx[i] ++;
+ }
+ }
+
+ void glp_set_multiobj_coef(glp_prob *lp, int objn, int j, double coef)
+ { glp_tree *tree = lp->tree;
+ int idx;
+ if (tree != NULL && tree->reason != 0)
+ xerror("glp_set_multiobj_coef: operation not allowed\n");
+ if (!(0 <= j && j <= lp->n))
+ xerror("glp_set_multiobj_coef: j = %d; column number out of range\n"
+ , j);
+ if (!(0< objn && objn <= lp->cobj_num+1))
+ xerror("glp_set_multiobj_coef: objn = %d; object number out of
range\n"
+ " use glp_set_multiobj_number() to set the number of
objectives\n"
+ , objn);
+ objn--;
+ if(objn == lp->cobj_num)
+ /* next object number, not seen yet */
+ { mk_multiobj_space(lp,objn); /* add space at the end */
+ lp->cobj_val[objn] = 0.0;
+ idx = lp->cobj_end;
+ lp->cobj_idx[objn] = idx;
+ mk_multiobj_space(lp,idx);
+ lp->cobj_idx[idx]=-1;
+ lp->cobj_num ++;
+ }
+ if (j == 0)
+ { lp->cobj_val[objn] = coef;
+ return;
+ }
+ for (idx=lp->cobj_idx[objn]; lp->cobj_idx[idx]>=0; idx++)
+ { if(lp->cobj_idx[idx] == j)
+ { lp->cobj_val[idx] = coef;
+ return;
+ }
+ }
+ xassert(lp->cobj_idx[idx]==-1);
+ if (coef == 0.0) /* do not store zero coeffs */
+ return;
+ mk_multiobj_space(lp,idx+1);
+ lp->cobj_idx[idx+1]=-1;
+ lp->cobj_val[idx+1] = lp->cobj_val[idx];
+ lp->cobj_idx[idx] = j;
+ lp->cobj_val[idx] = coef;
+ return;
+ }
+ #endif
+
/***********************************************************************
* NAME
*
***************
*** 1348,1353 ****
--- 1507,1516 ----
if (!(1 <= ncs && ncs <= lp->n))
xerror("glp_del_cols: ncs = %d; invalid number of columns\n",
ncs);
+ #ifdef CSL
+ if( lp->cobj_num > 0)
+ xerror("glp_del_cols: multiple objectives defined, cannot delete
columns\n");
+ #endif
for (k = 1; k <= ncs; k++)
{ /* take the number of column to be deleted */
j = num[k];
***************
*** 1451,1456 ****
--- 1614,1635 ----
dest->pbs_stat = prob->pbs_stat;
dest->dbs_stat = prob->dbs_stat;
dest->obj_val = prob->obj_val;
+ #ifdef CSL
+ if (prob->cobj_end > 0)
+ { dest->cobj_num = prob->cobj_num;
+ dest->cobj_max = prob->cobj_max;
+ dest->cobj_end = prob->cobj_end;
+ dest->cobj_idx = talloc(prob->cobj_max,int);
+ dest->cobj_val = talloc(prob->cobj_max,double);
+ memcpy(dest->cobj_idx,prob->cobj_idx,prob->cobj_end*sizeof(int));
+ memcpy(dest->cobj_val,prob->cobj_val,prob->cobj_end*sizeof(double));
+ }
+ else
+ { dest->cobj_end = 0;
+ dest->cobj_max = 0;
+ dest->cobj_num = 0;
+ }
+ #endif
dest->some = prob->some;
dest->ipt_stat = prob->ipt_stat;
dest->ipt_obj = prob->ipt_obj;
***************
*** 1564,1569 ****
--- 1743,1752 ----
if (lp->bfcp != NULL) xfree(lp->bfcp);
#endif
if (lp->bfd != NULL) bfd_delete_it(lp->bfd);
+ #ifdef CSL
+ if(lp->cobj_idx != NULL) xfree(lp->cobj_idx);
+ if(lp->cobj_val != NULL) xfree(lp->cobj_val);
+ #endif
return;
}
diff -r -c a/glpapi02.c b/glpapi02.c
*** a/glpapi02.c 2015-11-08 08:00:00.000000000 +0100
--- b/glpapi02.c 2015-12-02 14:33:39.536769618 +0100
***************
*** 388,393 ****
--- 388,429 ----
return j == 0 ? lp->c0 : lp->col[j]->coef;
}
+ #ifdef CSL
+ /***********************************************************************
+ * NAME
+ *
+ * glp_get_multiobj_coef - retrieve obj. coefficient or constant term
+ *
+ * SYNOPSIS
+ *
+ * double glp_get_obj_coef(glp_prob *lp, int objno, int j);
+ *
+ * RETURNS
+ *
+ * The routine glp_get_multiobj_coef returns the objective coefficient at
+ * j-th structural variable (column) of the specified problem object.
+ *
+ * If the parameter j is zero, the routine returns the constant term
+ * ("shift") of the objective function. */
+
+ double glp_get_multiobj_coef(glp_prob *lp, int objno, int j)
+ { int idx;
+ if (!(0 <= j && j <= lp->n))
+ xerror("glp_get_multiobj_coef: j = %d; column number out of range\n"
+ , j);
+ if (!(0 < objno && objno <= lp->cobj_num))
+ xerror("glp_get_multiobj_coef: objno = %d; object number out of
range\n"
+ , objno);
+ if (j == 0)
+ return lp->cobj_val[objno-1];
+ for (idx=lp->cobj_idx[objno-1]; lp->cobj_idx[idx]>=0; idx++)
+ { if (lp->cobj_idx[idx] == j)
+ return lp->cobj_val[idx];
+ }
+ return 0.0;
+ }
+ #endif
+
/***********************************************************************
* NAME
*
diff -r -c a/glpapi06.c b/glpapi06.c
*** a/glpapi06.c 2015-11-08 08:00:00.000000000 +0100
--- b/glpapi06.c 2015-12-03 16:13:56.155023543 +0100
***************
*** 498,503 ****
--- 498,506 ----
parm->out_frq = 500;
parm->out_dly = 0;
parm->presolve = GLP_OFF;
+ #ifdef CSL
+ parm->mobj = GLP_OFF;
+ #endif
return;
}
diff -r -c a/glpk.h b/glpk.h
*** a/glpk.h 2015-11-08 08:00:00.000000000 +0100
--- b/glpk.h 2015-12-03 16:12:46.743852966 +0100
***************
*** 134,139 ****
--- 134,142 ----
int out_frq; /* spx.out_frq */
int out_dly; /* spx.out_dly (milliseconds) */
int presolve; /* enable/disable using LP presolver */
+ #ifdef CSL
+ int mobj; /* flag enable/disable multiobjective */
+ #endif
double foo_bar[36]; /* (reserved) */
} glp_smcp;
***************
*** 299,304 ****
--- 302,312 ----
void glp_set_obj_dir(glp_prob *P, int dir);
/* set (change) optimization direction flag */
+ #ifdef CSL
+ void glp_set_multiobj_number(glp_prob *P, int nobjs);
+ /* set the number of extra objectives to the problem object */
+ #endif
+
int glp_add_rows(glp_prob *P, int nrs);
/* add new rows to problem object */
***************
*** 322,327 ****
--- 330,340 ----
void glp_set_obj_coef(glp_prob *P, int j, double coef);
/* set (change) obj. coefficient or constant term */
+ #ifdef CSL
+ void glp_set_multiobj_coef(glp_prob *P, int objno, int j, double coef);
+ /* set (change) additional object coefficient or constant term */
+ #endif
+
void glp_set_mat_row(glp_prob *P, int i, int len, const int ind[],
const double val[]);
/* set (replace) row of the constraint matrix */
***************
*** 397,402 ****
--- 410,420 ----
double glp_get_obj_coef(glp_prob *P, int j);
/* retrieve obj. coefficient or constant term */
+ #ifdef CSL
+ double glp_get_multiobj_coef(glp_prob *P, int objno, int j);
+ /* retrieve extra obj coefficient of constant term */
+ #endif
+
int glp_get_num_nz(glp_prob *P);
/* retrieve number of constraint coefficients */
Only in b/: libglpk.la
Only in b/: Makefile
diff -r -c a/prob.h b/prob.h
*** a/prob.h 2015-11-08 08:00:00.000000000 +0100
--- b/prob.h 2015-11-26 13:27:15.711006623 +0100
***************
*** 83,88 ****
--- 83,102 ----
AVL *c_tree;
/* column index to find columns by their names; NULL means this
index does not exist */
+ #ifdef CSL
+ int cobj_num;
+ /* actual number of extra objectives. Objectives are stored in an
+ index and value list. IDX[i] points to the start of the i'th
+ extra objective, VAL[i] is the shift (constant tag). Starting at
+ IDX[i] the column indices and objective values are stored; the
+ last item has value -1 */
+ int cobj_max, cobj_end;
+ /* allocated and used number of objective storage */
+ int *cobj_idx;
+ /* index of objective coefficients */
+ double *cobj_val;
+ /* value of the objectives */
+ #endif
/*--------------------------------------------------------------*/
/* basis factorization (LP) */
int valid;
diff -r -c a/simplex/spxprim.c b/simplex/spxprim.c
*** a/simplex/spxprim.c 2015-11-08 08:00:00.000000000 +0100
--- b/simplex/spxprim.c 2015-12-06 15:56:29.339492861 +0100
***************
*** 652,657 ****
--- 652,717 ----
skip: return;
}
+ #ifdef CSL
+ /***********************************************************************
+ * next_objective - set up the next objective
+ *
+ * After restriciting the free variables to the solution space,
+ * the next objective is computed into the array lp->c */
+
+ static void next_objective(glp_prob *P, struct csa *csa, int objno,
+ int map[/*1+m+n */])
+ { /* set up next next objective into lp->c */
+ SPXLP *lp = csa->lp;
+ int n = lp->n;
+ int m = lp->m;
+ double *c = lp->c;
+ double *d = csa->d;
+ double *l = lp->l;
+ double *u = lp->u;
+ int *head = lp->head;
+ char *flag = lp->flag;
+ double tol = csa->tol_dj, tol1 = csa->tol_dj1;
+ double eps,dir;
+ int objidx,j,k;
+ /* walk thu the list of non-basic variables */
+ for (j = 1; j <= n-m; j++)
+ { k = head[m+j]; /* x[k] = xN[j] */
+ if (l[k] == u[k]) /* xN[j] is fixed, skip it */
+ continue;
+ /* determine absolute tolerance eps[j] */
+ eps = tol + tol1 *(c[k] >= 0.0 ? +c[k] : -c[k]);
+ /* check if xN[j] will be fixed */
+ if (d[j] < -0.5*eps)
+ { /* xN[j] should be able to increase */
+ if (flag[j]) /* but its upper bound is active */
+ l[k] = u[k];
+ }
+ else if (d[j] > +0.5*eps)
+ { /* xN[j] should be able to decrease */
+ if (!flag[j] && l[k] != -DBL_MAX)
+ /* but its lower bound is active */
+ u[k] = l[k];
+ }
+ }
+ /* set up the new objective to lp->c */
+ dir = P->dir == GLP_MIN ? +1.0 : -1.0; /* get direction */
+ memset(c, 0, (n+1)*sizeof(double)); /* set c[] to zero */
+ for (objidx=P->cobj_idx[objno-1]; P->cobj_idx[objidx]>0; objidx++)
+ { j = P->cobj_idx[objidx];
+ k=map[m+j];
+ if (k == 0)
+ continue;
+ if (k < 0)
+ k = -k;
+ /* xassert(1 <= k && k <= lp->n); */
+ c[k] = dir * P->cobj_val[objidx] * P->col[j]->sjj;
+ }
+ /* and save the objective in csa->c */
+ memcpy(csa->c, c, (n+1)*sizeof(double));
+ }
+ #endif
+
/***********************************************************************
* spx_primal - driver to primal simplex method
*
***************
*** 992,997 ****
--- 1052,1060 ----
#endif
SPXSE se;
int ret, *map, *daeh;
+ #ifdef CSL
+ int objs;
+ #endif
/* build working LP and its initial basis */
memset(csa, 0, sizeof(struct csa));
csa->lp = &lp;
***************
*** 1075,1080 ****
--- 1138,1157 ----
csa->inv_cnt = 0;
/* try to solve working LP */
ret = primal_simplex(csa);
+ #ifdef CSL
+ /* go over further objectives; restrict the solution space and
+ compute the next objective */
+ if (parm->mobj == GLP_ON)
+ for (objs = 1; csa->p_stat == GLP_FEAS &&
+ csa->d_stat == GLP_FEAS && objs <= P->cobj_num; objs++)
+ { /* restrict search to the solution space and set up the new goal
*/
+ next_objective(P,csa,objs,map);
+ /* reduced costs are invalid */
+ csa->d_st = 0;
+ /* if (csa->se) csa->se->valid=0; */
+ ret = primal_simplex(csa);
+ }
+ #endif
/* return basis factorization back to problem object */
P->valid = csa->lp->valid;
P->bfd = csa->lp->bfd;
diff -r -c a/simplex/spydual.c b/simplex/spydual.c
*** a/simplex/spydual.c 2015-11-08 08:00:00.000000000 +0100
--- b/simplex/spydual.c 2015-12-06 15:56:29.339492861 +0100
***************
*** 174,181 ****
--- 174,183 ----
xassert(!flag[j]);
else if (l[k] == -DBL_MAX && u[k] != +DBL_MAX)
xassert(flag[j]);
+ #ifndef CSL /* this might not hold */
else if (l[k] == u[k])
xassert(!flag[j]);
+ #endif
}
return;
}
***************
*** 690,695 ****
--- 692,746 ----
skip: return;
}
+ #ifdef CSL
+ /***********************************************************************
+ * next_objective - set up the next objective
+ *
+ * After restriciting the free variables to the solution space,
+ * the next objective is computed into the array lp->c */
+
+ static void next_objective(glp_prob *P, struct csa *csa, int objno,
+ int map[/*1+m+n */])
+ { SPXLP *lp = csa->lp;
+ int n = lp->n;
+ int m = lp->m;
+ double *c = lp->c;
+ double *d = csa->d;
+ double *l = lp->l;
+ double *u = lp->u;
+ double tol = csa->tol_dj;
+ double tol1 = csa->tol_dj1;
+ int *head = lp->head;
+ double dir, eps;
+ int objidx,j,k;
+ /* change the bounds where the variable is fixed */
+ for (j = 1; j <= n-m; j++)
+ { k = head[m+j]; /* x[k] = xN[j] */
+ if (l[k] == u[k]) /* skip it */
+ continue;
+ eps = tol + tol1 * (c[k]>=0.0 ? +c[k] : -c[k]);
+ if (d[j] > +0.5*eps) /* lower bound is active */
+ { u[k] = l[k]; xassert(!lp->flag[j]); }
+ else if (d[j] < -0.5*eps) /* upper bound is active */
+ { l[k] = u[k]; xassert(lp->flag[j]); }
+ }
+ /* save lover and upper bounds */
+ memcpy(csa->l, l, (1+n) * sizeof(double));
+ memcpy(csa->u, u, (1+n) * sizeof(double));
+ /* set up the next objective to lp->c */
+ dir = P->dir == GLP_MIN ? +1.0 : -1.0; /* get direction */
+ memset(c,0,(n+1)*sizeof(double)); /* set c[] to zero */
+ for (objidx=P->cobj_idx[objno-1]; P->cobj_idx[objidx]>0; objidx++)
+ { j = P->cobj_idx[objidx];
+ k=map[m+j];
+ if (k == 0)
+ continue;
+ /* xassert(1 <= k && k <= lp->n); */
+ c[k] = dir * P->cobj_val[objidx] * P->col[j]->sjj;
+ }
+ }
+ #endif
+
/***********************************************************************
* spy_dual - driver to dual simplex method
*
***************
*** 1055,1060 ****
--- 1106,1114 ----
#endif
SPYSE se;
int ret, *map, *daeh;
+ #ifdef CSL
+ int objs;
+ #endif
/* build working LP and its initial basis */
memset(csa, 0, sizeof(struct csa));
csa->lp = &lp;
***************
*** 1154,1159 ****
--- 1208,1228 ----
csa->inv_cnt = 0;
/* try to solve working LP */
ret = dual_simplex(csa);
+ #ifdef CSL
+ /* go over further objectives; restrict the solution space and
+ compute the next objective */
+ if (parm->mobj == GLP_ON)
+ for (objs = 1; csa->p_stat == GLP_FEAS &&
+ csa->d_stat == GLP_FEAS && objs <= P->cobj_num; objs++)
+ { /* restrict search to the solution space */
+ next_objective(P,csa,objs,map);
+ csa->d_st = 0;
+ csa->phase = 0;
+ if(parm->pricing == GLP_PT_PSE)
+ csa->se->valid = 0;
+ ret = dual_simplex(csa);
+ }
+ #endif
/* return basis factorization back to problem object */
P->valid = csa->lp->valid;
P->bfd = csa->lp->bfd;
_______________________________________________
Help-glpk mailing list
[email protected]
https://lists.gnu.org/mailman/listinfo/help-glpk
