Trig Heuristic Method - Alpha V1.1 - C++ version - 11/18/15
Under The Artistic License 2.0, https://opensource.org/licenses/Artistic-2.0
This alpha version enables an increased number of contraints. If more than 100 constraints are required, please increase the array bound definition.
//Project starts here: special case when objective variable distance to the origin intersects the k-th constraints is not yet address here
//Two variable optimization project, which includes more than three constrains
#include <iostream> // std::cout, std::cin
#include <cmath> // std:: sqrt, std::pow, std::abs
using namespace std;
int i=0, j=0; // incremental operators
int a=0; // incremental operator in an array that is used to record the smallest angle for the shortest distance gamma, among equally small angles
int b=0; // incremental operator used to organize constraints bound by the k-th constraint; these constraints are below a threshold that requires them to considered in the set of primary constraints
int f=0; // incremental operators
int r=0; // operator set to sum the number of primary reference constraints, which excludes erased secondary set of reference constraints, and system constraints, which are established based on the k-th bounds
int q=0; // operator used to identify second reference constraint after the k-th one; this constraint is included in a primary set of constraint that are not bounded by the k-th constraint
int k=0; // operator used to identify first reference constraint, which yields the shortest distance gamma and smallest angle.
int kk=0; //resequenced k-th constraint
long double p=0; // operator used to subtract one array from another, and operator used towards computing the power function
int num=2; // number of variables in the objective function
int t=6; // number of constraints
int u=0; // increment operator used to assign value to array beyond the 6th referenced constraint
int g=0; // operator used to replace the primary set of reference constraints and erase secondary set of constraints
long double newconstr[100][2]; // matrix containing six of the original constraints and two project generated constraints
long double formatconstr[100][2]; // matrix containing set of primary reference constraints, which excludes erased secondary set of reference constraints
long double epi[2][2]; // matrix inverse of the matrix containing two constraints meant to be used for obtaining the quantity for each variable
long double x=0; // operator used to sum the squares, to measure distance between constraint coordinates and variable coefficients, by taking the square root for the sum of squares
int y=0; // identifier when the angle for the k-th constraint is null
int s=0; // identifier for goto loop, when angle for th k-th constraint is null
long double z=0; // operator used to sum the squares, to measure distance between variable coefficients and origin
int LC[100] = {};
long double LK[100] = {}; // measure distance between the k-th constraint and second likely constraint
long double L7[100] = {}; // array used to store angles
long double LA[100] = {}; // measure distance between constraint coordinates and origin; there are six constraints specified by the user and two system generated
// array used to store the distance between system constraint coordinates and origin; there are six constraints specified by the user and two system generated
long double L6[100] = {}; // measure distance between constraint coordinates and variable coefficients; there are six constraints specified by the user and two system generated
// array used to store the distance between system constraint coordinates and variable coefficients; there are six constraints specified by the user and two system generated
long double L8[100] = {}; // compute cosine, and thereafter record the angle, value of the triangle formed between the variable coefficients, a given constraint, and the origin
long double L9[100] = {}; // multiply distance between the variable coefficients and a given constraints and the cosine value; this provides distance between the two points if they were aligned between the variable coefficients and the origin
long double LB[100] = {}; // array that includes all constraints that will be excluded from the set of primary constraints, since the array includes all constraints bound by the k-th constraint below a target threshold
long double LE[100] = {}; // array used record the angle of the constraint yielding the smallest distance
int LF[100] = {}; // array used store which constraint yields the smallest distance
// array used to store the distance between system constraint coordinates and k-th constraint
long double w=0, v=0; // angles recorded from triangles formed, which are used to determine which plane formed by the binding a system constraint and the k-th constraint intersects the line formed by binding the coefficient variables and the origin
long double sol=0; // operator used to sum parameters for each constraint; this is equivalent to a sumproduct of the constraint with a unity vector
long double SOLL[2]; // array used to store the solution set for the maximization project
long double N[2] = { 18.3333333333333333, 20 }; // array of coefficients for variables in the objective function, which will be used in the maximization project within the prescribed constraint
// parameters for each of the six constraints are assigned to the variable to be optimized
long double constraints[][2]= {
{1, 2}, {3, 4}, {5, 6}, {7, 8}, {9, 10}, {11, 12}, {15, 12}, {12, 15},
{13, 16}, {14, 17}, {15, 18}, {16, 19}, {17, 20}, {18, 21}, {19, 22},
{20, 23}
};
// limiting parameters for each of the six constraints
long double limit[]= {
10, 20, 30, 40, 50, 60, 100, 100, 100, 100, 100, 100, 120, 120, 120, 120
};
//list of functions
void cosine_and_angle(long double array1[], long double array2[], int upper_bound); //compute cosine value and angle of the triangle formed between the variable coefficients, a given constraint, and the origin
void min(long double array3[], int upper_bound); //function used to determine the minimum within an array and below an upper bound
void matrix_inverse(long double array4[2][2], int amount6, int amoutn7); //function used to compute the inverse of a 2 x 2 matrix
long double cosine(long double amount1, long double amount2, long double amount3); //function used to compute the cosine
void dist_constr_var (long double array5[][2], int amount5, int amount4, long double return_array[], int upper_bound, int increment); //measure distance between constraint coordinates and variable coefficients, by taking the square root for the sum of squares
//int main begins
int main ()
{
// defined coefficients are provided for each variable of the objective function for the maximization problem
for(i=0;i<t;i++){
for(j=0;j<num;j++){
newconstr[i][j] = constraints[i][j]/ limit[i];
}
}
//measure distance between constraint coordinates and origin, by taking the square root for the sum of squares
dist_constr_var (newconstr, 0, 0, LA, t, 0) ;
//measure distance between constraint coordinates and variable coefficients, by taking the square root for the sum of squares
dist_constr_var (newconstr, 0, 1, L6, t, 0) ;
//measure distance between variable coefficients and origin, by taking the square root for the sum of squares
for(j=0;j<num;j++){
z += pow(N[j],2);
}
z = sqrt(z); //measure distance between variable coefficients and origin
//compute cosine value and angle of the triangle formed between the variable coefficients, a given constraint, and the origin
cosine_and_angle(LA,L6,t) ;
//find the smallest distance gamma recorded in array L9
min(L9, t);
//record all constraints for which the smallest distance gamma is equal
for(i=0;i<t;i++){
if (L9[i] == x) {
LF[a]=i; //array used record which constraint yields the smallest distance gamma
LE[a]=L8[i]; // array used record the angle of the constraint yielding the smallest distance gamma
a++;
}
}
//record the smallest angle for the shortest distance gamma, which is recorded in the set of ‘’a’’ for equally small angles
min(LE, a);
//record the k-th constraint that yields the shortest distance gamma and angle, which is the first reference constraint
for(i=0;i<a;i++){
if (LE[i] == x) {
k = LF[i];
x=0;
}
}
//determnine if the anlge for the k-th constraint is null
for(j=0;j<num;j++){
LE[j] = N[j]/ newconstr[k][j];
x += LE[j];
}
for(j=0;j<num;j++){
if ( x/num - LE[j] < 1e-10) { // y==1, this case applies when the angle for the k-th constraint is null, within rounding erro
p++;
}
}
if (p == num){
y = 1;
}
p = 0;
//mapping of back-up constraints by using the k-th constraint as binding
newconstr[6][0] = newconstr[k][0];
newconstr[6][1] = 0;
newconstr[7][0] = 0;
newconstr[7][1] = newconstr[k][1];
if (y != 1) {
//record the distance between the newly created system constraint, by taking the square root of the sum of squares
dist_constr_var (newconstr, 0, 0, LA, num, t) ;
dist_constr_var (newconstr, 0, 1, L6, num, t) ;
//measure distance between k-th constraint and system generated constraint
dist_constr_var (newconstr, 1, 0, LK, num, t);
//back-up constraints are now available
//measure the angle for each triangle and verify whether the angle formed by the k-th constraint is the difference between the former two
for(i=0;i<num;i++){
u = i+t;
v = cosine(LK[u], L6[u], L6[k]);
v = acos(v);
w = cosine(LA[u], z, L6[u]);
w = acos(w);
//determine which of the two system constraints should be used as the second reference L-th constraint with the k-th one
//computational error prevents the certitude of additional for angles; an error of 1e-17 is chosen as an acceptable threshold
if (abs( L8[k] + w - v )< 1e-10 ){
q=i+t;
}
}
for(j=0;j<num;j++){
if (newconstr[q][j] == newconstr[k][j]){
a = j;
v = newconstr[k][j];
}
}
}
kk=k;
repeat:
k=kk;
//verify whether all constraints, which exclude the k-th constraint, are bound by the k-th primary constraint and its shared axis with the q-th constraint
//repeat:
b=0;
g=0;
u = 2+t;
for(i=0;i<u;i++) {
f=0;
x=0;
// if and only if the angle for the k-th constraint is not null, then all constraint that are below a threshold, which is defined by the segment binding the k-th constraint and the q-th constraint, will chosen for removal
// review all constraints not excluded from the set of primary constraints, since the array includes all constraints bound by the k-th constraint below a target threshold
if (y != 1) {
if (i != k) {
if (newconstr[i][a] < v) {
f++;
}
if (newconstr[i][a] >= v) {
x++;
p++;
}
}
if (i == k) {
x++;
}
}
if (y == 1) {
if (i != k) {
if (newconstr[i][s] < newconstr[k][s]) {//repeat loop for each "s" coordinate of the k-th constraint
f++;
}
if (newconstr[i][s] >= newconstr[k][s]) {//repeat loop for each "s" coordinate of the k-th constraint
x++;
p++;
}
}
if (i == k) {
x++;
}
}
if ( f != 0 ) {
LB[b] = i; //i-th constraint will be recorded in the LB array, which includes all constraints that will be excluded from the set of primary constraints
b++;// ”b” operator indicates the total number of constraints bound by the k-th constraint
}
if ( x != 0 ) {// array used to sort the reference set of primary constraints
LF[g] = i;
if (i == k) {
k=g;
}
g++;
}
}
f=0;
// if the number of constraints excluded from the set of primary constraints amounts to the total number constraints minus one, which is the k-th constraint, then the k-th constraint is the only binding one
a=0;
g=0;
if (p == 0) { a=1; } // a is an operator used to determine whether the first reference constraint is the only binding constraint
p=0;
if (a == 1) {
//first constraint reference provides bounds for solution; the L-th and k-th constraints are used to determine the solution
//the following provide the inverse of a 2x2 matrix and output of solution
if (y != 1) {
matrix_inverse(newconstr, q, k);
}
if (y == 1) {
for(i=0;i<num;i++){
q=i+t;
matrix_inverse(newconstr, q, k);
cout<< "\n";
}
}
s++;
if (s < num) {
goto repeat;
}
return 0; //the project ends here, if a=1
}
//matrix containing set of primary reference constraints, which excludes erased secondary set of reference constraints
r = (u-b);
for (i=0;i<r;i++){
a=LF[i];
LE[i]=LF[i];
for (j=0;j<num;j++){
formatconstr[i][j] = newconstr[a][j];
}
}
//measure distance between constraint coordinates and origin, by taking the square root for the sum of squares
dist_constr_var (formatconstr, 0, 0, LA, r, 0) ;
//measure distance between constraint coordinates and variable coefficients, by taking the square root for the sum of squares
dist_constr_var (formatconstr, 0, 1, L6, r, 0) ;
//recalculating cosine values
cosine_and_angle(LA,L6,r) ;
//begin to locate second reference constraint
dist_constr_var (formatconstr, 1, 0, LK, r, 0);
f=0;
for(i=0;i<r;i++){
if (i != k) {
// angle recorded from the reference k-th constraint and likely second constraint
v = cosine(LK[i], L6[k], L6[i]);
v = acos(v); // determine angle v
//sum of angles should amount to the addition of the k-th constraint angle and likely second reference constraint;
//this provides an array of distances between the k-th constraints and other constraints; these distances must intersect with the distance from the origin to the objective variable coordinates
//computational error prevents the certitude of additional for angles; an error of 1e-12 is chosen as an acceptable threshold rather than 1e-17
// L6 array used to store the second reference constraint after the k-th one; this constraint is included in a primary set of constraints that are not bounded by the k-th constraint
if (abs( L8[k] + L8[i] - v ) < 1e-10) {
L6[f] = i;
f++;
}
}
}
//determine angle formed between a slope orthogonal to gamma and the distance between the k-th constraint and the likely second reference constraint
for(i=0;i<f;i++){
a = L6[i];
L7[i] = acos( cosine( L9[a], L9[k], LK[a] ) );
}
//determine the minimum
min(L7, f);
for(i=0;i<f;i++){
if (L7[i]==x){
q=L6[i];
a=q;
}
}
if (y != 1) {
matrix_inverse(formatconstr, q, k); //the following provide the inverse of a 2x2 matrix and output solution
x=0;
p=0;
for(i=0;i<r;i++) {
if (i != q) {
if (i != k) {
for (j=0;j<num;j++) {
x += formatconstr[i][j] * SOLL[j];
if (1 - x < 1e-10) {
LF[i]=i;
p++;
}
}
x=0;
}
}
}
for(g=0;g<p;g++) {
q=LF[g];
matrix_inverse(formatconstr, q, k);
}
}
if (y == 1) {//case when the angle for the k-th constraint is null
matrix_inverse(formatconstr, q, k); //the following provide the inverse of a 2x2 matrix and output solution
x=0;
p=0;
for(i=0;i<r;i++) {
if (i != q) {
if (i != k) {
for (j=0;j<num;j++) {
x += formatconstr[i][j] * SOLL[j];
if (1 - x < 1e-10) {
LF[i]=i;
p++;
}
}
x=0;
}
}
}
for(g=0;g<p;g++) {
q=LF[g];
matrix_inverse(formatconstr, q, k);
}
s++;
if (s < num) {
goto repeat;
}
}
return 0;
//This is the end of the project
}
void cosine_and_angle(long double array1[], long double array2[], int upper_bound)
{
for(i=0; i < upper_bound ;i++){
//compute cosine value of the triangle formed between the variable coefficients, a given constraint, and the origin
L8[i] = cosine(array1[i], array2[i], z);
//multiply distance, which is measured between the variable coefficients and a given constraints, and the cosine value; this provides distance between the two points if they were aligned between the variable coefficients and the origin; this distance will be referenced as gamma
L9[i]=L8[i]*array2[i];
//the angle for the cosine
L8[i] = acos(L8[i]);
}
}
void min(long double array3[], int upper_bound)
{
//x will be assigned the minimum value
x = array3[0];
for(i=0 ; i < upper_bound ; i++) {
if (array3[i] < x) {
x = array3[i];
}
}
}
void matrix_inverse(long double array4[2][2], int amount6, int amount7) // inverse of 2 x 2 matrix and output of solution
{
x=0;
b=amount6;
a=amount7;
// matrix containing two constraints meant to be used for obtaining the quantity for each variable
//matrix inverse
//determinant
v = array4[a][0] * array4[b][1] - array4[a][1] * array4[b][0];
epi[0][0] = array4[b][1] / v;
epi[1][1] = array4[a][0] / v;
epi[0][1] = - array4[a][1] / v;
epi[1][0] = - array4[b][0] / v;
//output of solution
for(i=0;i<num;i++){
for(j=0;j<num;j++){
x += epi[i][j];
}
SOLL[i]=x;
x=0;
}
sol=0;
for(i=0;i<num;i++){
sol += N[i] * SOLL[i];
}
//optimal quantities for each variable are
cout<< "\noptimal quantities for each decision variable are respectively "<<SOLL[0]<<" "<< SOLL[1]<< "\n";
//optimal solution is
cout<< "optimal solution is " <<sol << "\n";
//set of optimal constraints which can be combined to determine the optimal quantity for decision variables
cout<< "\nthis set of constraints yields the optimal quantity for the decision variables \n";
b = LE[b]; //warning: convert double to int
a = LE[a]; //warning: convert double to int
cout<< "this k-th constraint is pre-requisit for any combination with the set \n[";
for(j=0;j<num;j++){
cout<<constraints[a][j]<<" ";
}
cout<<" <= "<<limit[a]<<"\n\n";
if ( b > t-1 ) {
cout<< "all other constraints within the set are listed \n";
for(j=0;j<num;j++){
cout<<newconstr[b][j]*limit[a]<<" ";
}
cout<<"] <= "<<limit[a]<<"\n\n";
}
if ( b <= t-1 ) {
cout<< "all other constraints within the set are listed \n";
for(j=0;j<num;j++){
cout<<constraints[b][j]<<" ";
}
cout<<"] <= "<<limit[b]<<"\n\n";
}
}
long double cosine(long double amount1, long double amount2, long double amount3) //function used to compute the cosine
{
return (pow(amount1,2) - pow(amount2,2) - pow(amount3,2)) / (-2 * amount2 * amount3);
}
//measure distance between constraint coordinates and variable coefficients, by taking the square root for the sum of squares
void dist_constr_var (long double array5[][2], int amount5, int amount4, long double return_array[], int upper_bound, int increment)
{
x=0;
for(i=0; i < upper_bound ;i++){
u = i + increment;
for(j=0;j<num;j++){
//measure distance between constraint coordinates and variable coefficients, by taking the square root for the sum of squares
x += pow((array5[u][j] - amount5*array5[k][j] - amount4*N[j]), 2);
}
// measure distance between constraint coordinates and variable coefficients
return_array[u] = sqrt(x);
x = 0;
}
}