public static void OptimByConvexDimDecomposition(string OptimProblemTag, string MaxOrMin, double[] C, double[,] A, string[] ComparatorVctr, double[] b, ref double[] x, ref string OptimalValue, ref string DominantDimension, ref string ActiveKnstraints, ref string RedundantKnstraints, ref string UpperUnBoundDim, ref string LowerUnBoundDim, ref double TimeToSolve_mSecs, ref string Comments) { // OptimByConvexDimDecomposition: by Gerardo L. Febres. www.gfebres.net. Version: 2019.07.29 ------------------------------------------ // // Solves a linear (max) optimization problem by recognizing the active constraints // // C A double[] ObjFunCoefsVectSTRC, double[,] ConstrCoefMatrxSTRC, // Comparator string[] ConstrComparatorVectSTRC, // b x double[] ResourcesVectSTRC, ref double[] OptimalCoordinates // // Returns the optimal point for a linear optimization problem using the Febres Algorith"m: // // x List of Optimal Coordinate values separated by the symbols "]0[" // g Dominant dimension. // ActiveKnstraints Structure with Active Constraints indexes separated by the symbols "]0[" and "]1[" as follows: // Dim 0 ]1[ ActvUpprConstr j ]0[ Dim 1 ]1[ ActvUpprConstr k ]0[ ...Dim n ]1[ ActvUpprConstr m // // Non-UpperBound Dimensions Structure with Non-UpperBound Dimension indexes separated by the symbols "]0[" as follows: // Dim 0 ]0[ Dim j ]0[ Dim k ]0[ ...]0[ Dim p // // Non-LowerBound Dimensions Structure with Non-LowerBound Dimension indexes separated by the symbols "]0[" as follows: // Dim 0 ]0[ Dim j ]0[ Dim k ]0[ ...]0[ Dim p // // RedundantKnstraints List of Redundant Constraints indexes separated by the symbols "]0[" // TimeToSolve_mSecs Time to solve in 0.00001 sec // Comments Comments // ------------------------------------------------------------------------------------------------------------------------ // // Explicit C sharp code for the Linear Optimization algorithm based on active constraints recognition. // The algorithm is not iterative and its complexity is comparable to the Simplex algorithm. // // To have the code working, copying it and pasting on Visual Studio may work. // // Search "A non-iterative method for optimizing linear functions in convex linearly constrained spaces" // by Gerardo L.Febres for the theoretical basis of this algorithm. // // ------------------------------------------------------------------------------------------------------------------------ // 3.1. Define and dimension set variables and arrays. --------------------------------------------------------------- // 3.1.1 Check dimensionality coherence. Define Variables and arrays. ------------------------------------------------- bool DimensionalityTestsPassed = true; int[] UpperUnBoundDimIndx = new int[0]; int[] LowerUnBoundDimIndx = new int[0]; if (MaxOrMin != null && MaxOrMin != "") { if (MaxOrMin.ToUpper().StartsWith("MAX")) MaxOrMin = "max"; else if (MaxOrMin.ToUpper().StartsWith("MIN")) MaxOrMin = "min"; else DimensionalityTestsPassed = false; } else { MaxOrMin = "max"; } if (C.GetUpperBound(0) != A.GetUpperBound(1)) { MessageBox.Show("Coeficient Matrix A and Objective Function Vector C must have the same number of columns." + "\n\r\n\r" + "Please check Problem:" + "\n\r\n\r" + OptimProblemTag, "MsgCode: 2017.08.17.A", MessageBoxButtons.OK, MessageBoxIcon.Information); DimensionalityTestsPassed = false; } if (b.GetUpperBound(0) != A.GetUpperBound(0)) { MessageBox.Show("Coeficient Matrix A and Resource Vector b must have the same number of elements." + "\n\r\n\r" + "Please check Problem:" + "\n\r\n\r" + OptimProblemTag, "MsgCode: 2017.08.17.B", MessageBoxButtons.OK, MessageBoxIcon.Information); DimensionalityTestsPassed = false; } if (ComparatorVctr.GetUpperBound(0) != A.GetUpperBound(0)) { MessageBox.Show("Coeficient Matrix A and Vector ComparatorVctr must have the same number of elements." + "\n\r\n\r" + "Please check Problem:" + "\n\r\n\r" + OptimProblemTag, "MsgCode: 2017.08.17.C", MessageBoxButtons.OK, MessageBoxIcon.Information); DimensionalityTestsPassed = false; } if (b.GetUpperBound(0) < A.GetUpperBound(1)) { MessageBox.Show("Any PPL in this resentation must have at least al many Constraints as Dimensions." + "\n\r\n\r" + "Please check Problem: " + "\n\r\n\r" + OptimProblemTag, "MsgCode: 2017.10.19.A", MessageBoxButtons.OK, MessageBoxIcon.Information); DimensionalityTestsPassed = false; } // Solve if (DimensionalityTestsPassed) { TimeToSolve_mSecs = Convert.ToDouble(DateTime.Now.ToString("ss.fffffff")); // Define Variables and arrays. double BigN = 999999999; double Tiny = 0.00000001; int m = A.GetUpperBound(0) + 1; // Number of Constraints m int n = A.GetUpperBound(1) + 1; // Number of dimensions n // Organize Inward Constraits first and then outward. if (MaxOrMin == "max") { } else { // (MaxOrMin == "min") MessageBox.Show("Code is to be extended for minimization problems '" + "\n\r\n\r" + "Try convertting problem '" + OptimProblemTag + "' into a maximization problem.", "MsgCode: 2019.04.07.C", MessageBoxButtons.OK, MessageBoxIcon.Information); DimensionalityTestsPassed = false; } // Reference Vectors Vof and Vctr double[] Vof = new double[n]; // Obj.Funct.Unitary Director Vector double VofNorm = 0; for (int i = 0; i <= n - 1; i++) { Vof[i] = C[i]; VofNorm = VofNorm + Math.Pow(Vof[i], 2); } VofNorm = Math.Pow(VofNorm, 0.5); for (int i = 0; i <= n - 1; i++) { Vof[i] = Vof[i] / VofNorm; } // Constraints' Director Vectors considering type of constraints. double[,] Aji = new double[A.GetUpperBound(0) + 1, A.GetUpperBound(1) + 1]; for (int j = 0; j <= Aji.GetUpperBound(0); j++) { for (int i = 0; i <= Aji.GetUpperBound(1); i++) Aji[j, i] = A[j, i]; } string[] ConstraintType = new string[m]; double[,] Vji = new double[m, n]; for (int j = 0; j <= m - 1; j++) { if (ComparatorVctr[j] == "<=" || ComparatorVctr[j] == "<") ConstraintType[j] = "Inward"; else { if (ComparatorVctr[j] == ">=" || ComparatorVctr[j] == ">") { ConstraintType[j] = "Outward"; for (int i = 0; i <= n - 1; i++) Aji[j, i] = -Aji[j, i]; } else { if (ComparatorVctr[j] == "=") ConstraintType[j] = "Equality"; else { MessageBox.Show("The comparator '" + ComparatorVctr[j] + "' cannot be processed." + "\n\r\n\r" + "Check the comparators for: " + "\n\r\n\r" + OptimProblemTag, "MsgCode: 2018.12.22.A", MessageBoxButtons.OK, MessageBoxIcon.Information); DimensionalityTestsPassed = false; } } } double VjNorm = 0; for (int i = 0; i <= n - 1; i++) { if(A[j, i] == 0) { Vji[j, i] = 0; } else { Vji[j, i] = 1 / A[j, i]; } VjNorm = VjNorm + Math.Pow(Vji[j, i], 2); } VjNorm = Math.Pow(VjNorm, 0.5); for (int i = 0; i <= n - 1; i++) Vji[j, i] = Vji[j, i] / VjNorm; } if (DimensionalityTestsPassed) { // 3.1.2. Determine dominant dimension g. Fastest Growing Dimension g. -------------------------------------------- int g = 0; double MaxC = -BigN; for (int j = 0; j <= m - 1; j++) { if (ConstraintType[j] == "Inward") { for (int i = 0; i <= n - 1; i++) { if (A[j, i] != 0) { if (C[i] / A[j, i] > MaxC) { g = i; MaxC = C[i] / A[j, i]; } } } } } DominantDimension = g.ToString(); // Constraint-Objective-Function-Resource Impact Index. Constraint Coeff. Parameter E[j] (3.1.2.a)------------------- double[] E = new double[m]; for (int j = 0; j <= m - 1; j++) { E[j] = A[j, 0]; for (int i = 1; i <= n - 1; i++) E[j] = E[j] + A[j, i] * C[i] / C[0]; } // Constraint Closest Points By the objective-function direction d[j]. (3.1.2.b)-------------------------------------- double[] d = new double[m]; // d[j] Dist. from origin to Constraint j by obj.-funct. direction (BOFD) for (int j = 0; j <= m - 1; j++) { d[j] = 1; for (int i = 1; i <= n - 1; i++) { d[j] = d[j] + Math.Pow((C[i] / C[0]), 2); } d[j] = b[j] / E[j] * Math.Pow(d[j], 0.5); if(ConstraintType[j] == "Outward" && d[j] < 0) {d[j] = -d[j];} } // 3.1.3. Determine the coordinates Xcp for each dimension i of closest point for all Constraints j. Xcp[j, i]-------- double[,] Xcp = new double[m, n]; // Xcp[j, i] Coordinates of Constraint's closest point by obj.-funct. direction for (int j = 0; j <= m - 1; j++) { Xcp[j, 0] = b[j] / E[j]; if (ConstraintType[j] == "Outward" && Xcp[j, 0] < 0) { Xcp[j, 0] = -Xcp[j, 0]; } for (int i = 1; i <= A.GetUpperBound(1); i++) Xcp[j, i] = C[i] / C[0] * Xcp[j, 0]; } // 3.1.4. Determine the Objective-Function's Angles and the On Principal Planes Projected Constraint Angles. ------------------------------------------------------ double[] APjctnOi = new double[n]; // APjctnOi[,] angle between Projection of Obj.Funct. vector on Plane i-g and axis g double[,] APjctGji = new double[m, n]; // APjctGji[,] angle between Projection of Constraint j's Normal vector on Plane i-g and axis g for (int j = 0; j <= m - 1; j++) { for (int i = 0; i <= n - 1; i++) { if (i != g) { APjctnOi[i] = Math.Atan(C[i] / C[g]); if (Aji[j, i] == 0 && Aji[j, g] == 0) { // Constraint j is paralell to plane ig APjctGji[j, i] = BigN; } else if (Aji[j, i] >= 0 && Aji[j, g] >= 0) { APjctGji[j, i] = Math.Atan(Aji[j, i] / Aji[j, g]); if (Aji[j, g] == 0) { APjctGji[j, i] = Math.Sign(Aji[j, i]) * APjctGji[j, i]; } } else if (Aji[j, i] >= 0 && Aji[j, g] < 0) { APjctGji[j, i] = Math.Atan(Aji[j, i] / Aji[j, g]); if (Aji[j, i] == 0) { APjctGji[j, i] = - Math.PI + APjctGji[j, i]; } else { APjctGji[j, i] = Math.PI + APjctGji[j, i]; } } else if (Aji[j, i] < 0 && Aji[j, g] >= 0) { APjctGji[j, i] = Math.Atan(Aji[j, i] / Aji[j, g]); if (Aji[j, g] == 0) { APjctGji[j, i] = -Math.PI / 2; } } else if (Aji[j, i] < 0 && Aji[j, g] < 0) { APjctGji[j, i] = Math.Atan(Aji[j, i] / Aji[j, g]) - Math.PI ; } } } } // 3.2. Study problem's upper constraints.--------------------------------------------------------------------- // 3.2.1. Recognize Upper-Constraints for all projections over the principal planes. Apply Equation (3).-------------- bool[,] UpperBoundsDIMg = new bool[m, n]; for (int j = 0; j <= m - 1; j++) { for (int i = 0; i <= n - 1; i++) { if(i != g) { double SignedAlphPrjctDif = APjctGji[j, i] - APjctnOi[i]; if (SignedAlphPrjctDif > Math.PI / 2) { } // DoesNotUpperBindDIMg[j, i] = true;} else if (SignedAlphPrjctDif >= 0) {UpperBoundsDIMg[j, i] = true;} else if (SignedAlphPrjctDif < -Math.PI / 2) { } // DoesNotUpperBindDIMg[j, i] = true; } } } // Review upper-Constraints'Bound-conditions. bool ProblemBoundnessCOndition = true; for (int i = 0; i <= n - 1; i++) { if (i != g) { int j = 0; string ActiveUpperKnstraints = ""; for (j = 0; j <= m - 1; j++) { if (UpperBoundsDIMg[j, i]) { if (ActiveUpperKnstraints == "") ActiveUpperKnstraints = i + "]1[" + j; else ActiveUpperKnstraints = ActiveUpperKnstraints + "]0[" + i + "]1[" + j; break; } } if (j > m - 1) { ProblemBoundnessCOndition = false; Array.Resize(ref UpperUnBoundDimIndx, UpperUnBoundDimIndx.GetUpperBound(0) + 2); UpperUnBoundDimIndx[UpperUnBoundDimIndx.GetUpperBound(0)] = i; } } } if (ProblemBoundnessCOndition) { // 3.2.2 Identify the upper-constraint with the shortest hypotenuse for each principal plane projection. Apply Equation (7). double BigPostv = BigN; double[,] uProjFactr = new double[m, n]; // uProjFactr[j, i] Rel. of cuts to axes i and g of prjectn. on princ. plane i-g for Constraint j int[] uShrtstHypo = new int[n]; // uShrtstHypo[j, i] index of Constraint j with Shortest Hipo in principal plane i-g. double[] uShrtstHypoSize = new double[n]; // uShrtstHypoSize[j, i] Size of Constraint j with Shortest Hipo in principal plane i-g. for (int i = 0; i <= n - 1; i++) uShrtstHypo[i] = -1; for (int i = 0; i <= n - 1; i++) { if (i != g) { uShrtstHypoSize[i] = BigPostv + 1; for (int j = 0; j <= m - 1; j++) { if (UpperBoundsDIMg[j, i]) { // 3.2.2.a Compute the Projection-Factor for upper-constraints. Equation (5).------------------------ int dimi = 0; for (dimi = 0; dimi <= n - 1; dimi++) { if (dimi != i && dimi != g) { if (A[j, dimi] != 0) { break; } } } if (dimi > n - 1) { uProjFactr[j, i] = 1; } else { uProjFactr[j, i] = (Xcp[j, g] * A[j, g] + Xcp[j, i] * A[j, i]) / b[j]; if (uProjFactr[j, i] < 0) { uProjFactr[j, i] = -uProjFactr[j, i]; } } // ------------------------------------------------------ double uHypoSize = BigPostv; if (A[j, i] == 0 || A[j, g] == 0) { if (uShrtstHypo[i] < 0) { // Select Constraint if is the first uShrtstHypo[i] = j; uShrtstHypoSize[i] = uHypoSize; } else { // Hypo is Infinite..Make Decision if (A[j, g] == 0) { if (uShrtstHypoSize[i] == BigPostv) { // Select the containing closest to the origin double IuCut = uProjFactr[j, i] * b[j] / A[j, i]; double IuhCut = uProjFactr[j, i] * b[uShrtstHypo[i]] / A[uShrtstHypo[i], i]; if (IuCut < IuhCut) { uShrtstHypo[i] = j; uShrtstHypoSize[i] = uHypoSize; } } else { // Keep current uShrtstHypoSize[i] } } } } else { double IuCut = uProjFactr[j, i] * b[j] / A[j, i]; double GuCut = uProjFactr[j, i] * b[j] / A[j, g]; uHypoSize = Math.Pow(Math.Pow(IuCut, 2) + Math.Pow(GuCut, 2), 0.5); if (uHypoSize < uShrtstHypoSize[i]) { uShrtstHypo[i] = j; uShrtstHypoSize[i] = uHypoSize; } else { if (Math.Abs(uHypoSize - uShrtstHypoSize[i]) < Tiny) { // Constrs. j and uShrtstHypoSize[i] are equally sized. MessageBox.Show("It does enter here!!! 2019.03.29.A"); // Select the one closer in the dimensions orthogonal to axis i and g if (uHypoSize == BigPostv) { if (GuCut < uShrtstHypoSize[i]) { uShrtstHypo[i] = j; uShrtstHypoSize[i] = uHypoSize; } else { double TieMeasrShrtsJ = 0; double TieMeasrThis = 0; for (int TieUIdx = 0; TieUIdx <= n - 1; TieUIdx++) { TieMeasrThis = TieMeasrThis + b[j] / A[j, TieUIdx]; TieMeasrShrtsJ = TieMeasrShrtsJ + b[uShrtstHypo[i]] / A[uShrtstHypo[i], TieUIdx]; } if (TieMeasrThis < TieMeasrShrtsJ) { uShrtstHypo[i] = j; uShrtstHypoSize[i] = uHypoSize; } else { // Two Constraints are overlapping } } } else { // Identical intersection with plane i-g. Take the closer to the origin if (Math.Abs(d[i] - d[uShrtstHypo[i]]) < Tiny) { // Identical Constraints. Ignore the most recent } else { if (d[j] < d[uShrtstHypo[i]]) { uShrtstHypo[i] = j; uShrtstHypoSize[i] = uHypoSize; } } } } } } } } } } // 3.2.3 Build array with the possible active upper-constraints in each projections over the principal planes. // Rank Upper constraints according to the Constraint Gradient criterion. ---------------------------- int[][] RankedUpperConstrJndx = new int[n][]; double[][] CutOverAxisI = new double[n][]; double[][] CutOverAxisG = new double[n][]; double[][] UpperConstrGrad = new double[n][]; string[][] OverlapCnstrSetU = new string[n][]; for (int i = 0; i <= n - 1; i++) { if (i != g) { Array.Resize(ref RankedUpperConstrJndx[i], 0); Array.Resize(ref CutOverAxisI[i], 0); Array.Resize(ref CutOverAxisG[i], 0); Array.Resize(ref UpperConstrGrad[i], 0); Array.Resize(ref OverlapCnstrSetU[i], 0); for (int j = 0; j <= m - 1; j++) { if (UpperBoundsDIMg[j, i]) { double Iuh = uProjFactr[uShrtstHypo[i], i] * b[uShrtstHypo[i]] / A[uShrtstHypo[i], i]; double Guh = uProjFactr[uShrtstHypo[i], i] * b[uShrtstHypo[i]] / A[uShrtstHypo[i], g]; bool includeConstrUj = false; double Gradig = 0; if (j == uShrtstHypo[i]) { // this is the Minimum Hypothenuse constraint. Include. includeConstrUj = true; Gradig = A[j, g] / A[j, i]; } else { if (A[j, i] == 0 || A[j, g] == 0) { // The plane j is parallel or perpendicular to plane ig and hypo infinite if (A[j, i] == 0) { // Plane j is parallel to axis i. Hypo is infinite. double Gu = uProjFactr[j, i] * b[j] / A[j, g]; if (Gu < Guh) { includeConstrUj = true; Gradig = BigPostv; } } if (A[j, g] == 0) { // Plane j is parallel to axis g. Hypo is infinite. double Iu = uProjFactr[j, i] * b[j] / A[j, i]; if (Iu < Iuh) { includeConstrUj = true; Gradig = 0; } } } else { double Iu = uProjFactr[j, i] * b[j] / A[j, i]; double Gu = uProjFactr[j, i] * b[j] / A[j, g]; double IntersectDetect = (Iu - Iuh) * (Gu - Guh); // if (IntersectDetect == 0) { // Degenerative intersection. Select the one with shortest hypotenuse. if (IntersectDetect >= -0.0005 && IntersectDetect <= 0.0005) { // Degenerative intersection. Select the one with shortest hypotenuse. includeConstrUj = true; Gradig = A[j, g] / A[j, i]; } else { if (IntersectDetect < 0) { // Constraint j cuts constraint uShrtstHypo[i] includeConstrUj = true; Gradig = A[j, g] / A[j, i]; } } } } if (includeConstrUj) { double AxisICutValue = b[j] / A[j, i]; double AxisGCutValue = b[j] / A[j, g]; int DegeneraIdx = 0; for (DegeneraIdx = 0; DegeneraIdx <= RankedUpperConstrJndx[i].GetUpperBound(0); DegeneraIdx++) { if(AxisICutValue == CutOverAxisI[i][DegeneraIdx] && AxisGCutValue == CutOverAxisG[i][DegeneraIdx]) { // Including this Constraint will cause Deneration of the Equation System. 2019.07.20 break; } } if(DegeneraIdx <= RankedUpperConstrJndx[i].GetUpperBound(0) || RankedUpperConstrJndx[i].GetUpperBound(0) < 0){ if (RankedUpperConstrJndx[i].GetUpperBound(0) < 0){ Array.Resize(ref OverlapCnstrSetU[i], OverlapCnstrSetU[i].GetUpperBound(0) + 2); OverlapCnstrSetU[i][OverlapCnstrSetU[i].GetUpperBound(0)] = j.ToString(); } else{ OverlapCnstrSetU[i][DegeneraIdx] = OverlapCnstrSetU[i][DegeneraIdx] + ","+ j.ToString(); } } else{ Array.Resize(ref OverlapCnstrSetU[i], OverlapCnstrSetU[i].GetUpperBound(0) + 2); OverlapCnstrSetU[i][OverlapCnstrSetU[i].GetUpperBound(0)] = j.ToString(); } Array.Resize(ref RankedUpperConstrJndx[i], RankedUpperConstrJndx[i].GetUpperBound(0) + 2); Array.Resize(ref CutOverAxisI[i], CutOverAxisI[i].GetUpperBound(0) + 2); Array.Resize(ref CutOverAxisG[i], CutOverAxisG[i].GetUpperBound(0) + 2); Array.Resize(ref UpperConstrGrad[i], UpperConstrGrad[i].GetUpperBound(0) + 2); int AtJuIdx = 0; while (Gradig < UpperConstrGrad[i][AtJuIdx] && AtJuIdx < RankedUpperConstrJndx[i].GetUpperBound(0)) { AtJuIdx++; } for (int SldIdx = RankedUpperConstrJndx[i].GetUpperBound(0); SldIdx > AtJuIdx; SldIdx--) { RankedUpperConstrJndx[i][SldIdx] = RankedUpperConstrJndx[i][SldIdx - 1]; CutOverAxisI[i][SldIdx] = CutOverAxisI[i][SldIdx - 1]; CutOverAxisG[i][SldIdx] = CutOverAxisG[i][SldIdx - 1]; UpperConstrGrad[i][SldIdx] = UpperConstrGrad[i][SldIdx - 1]; } RankedUpperConstrJndx[i][AtJuIdx] = j; CutOverAxisI[i][AtJuIdx] = AxisICutValue; CutOverAxisG[i][AtJuIdx] = AxisGCutValue; UpperConstrGrad[i][AtJuIdx] = Gradig; } } } } } // 3.3 Study the 2D-projected problem’s lower-constraints.--------------------------------------------------- // 3.3.1 Recognize Lower-Constraints for all projections over the principal planes. Apply Equation (4).-------- bool[,] LowerBoundsDIMg = new bool[m, n]; for (int j = 0; j <= m - 1; j++) { for (int i = 0; i <= n - 1; i++) { if (i != g) { for (int u = 0; u <= 0; u++) { // or (int u = 0; u <= RankedUpperConstrJndx[i].GetUpperBound(0); u++) if (APjctGji[j, i] <= -Math.PI) { APjctGji[j, i] = APjctGji[j, i] + 2 * Math.PI; } if (APjctGji[RankedUpperConstrJndx[i][u], i] <= -Math.PI) { APjctGji[RankedUpperConstrJndx[i][u], i] = APjctGji[RankedUpperConstrJndx[i][u], i] + 2 * Math.PI; } double SignedAlphPrjctDif = APjctGji[j, i] - APjctGji[RankedUpperConstrJndx[i][u], i]; if (SignedAlphPrjctDif > Math.PI) { SignedAlphPrjctDif = 2 * Math.PI - SignedAlphPrjctDif; } if (SignedAlphPrjctDif < - Math.PI) { SignedAlphPrjctDif = SignedAlphPrjctDif + 2 * Math.PI; } if (SignedAlphPrjctDif >= 0) { } // DoesNotLowerBindDIMg[j, i] = true; else if (SignedAlphPrjctDif <= -Math.PI) { } // DoesNotLowerBindDIMg[j, i] = true; else { LowerBoundsDIMg[j, i] = true; } } } } } for (int i= 0; i <= n - 1; i++) { // Review Lower-Constraints'Bound-conditions if (i != g) { int j = 0; string ActiveLowerKnstraints = ""; for (j = 0; j <= m - 1; j++) { if (LowerBoundsDIMg[j, i]) { if(ActiveLowerKnstraints=="") ActiveLowerKnstraints = i + "]1[" + j; else ActiveLowerKnstraints = ActiveLowerKnstraints + "]0[" + i + "]1[" + j; break; } } if (j > m - 1) { ProblemBoundnessCOndition = false; Array.Resize(ref LowerUnBoundDimIndx, LowerUnBoundDimIndx.GetUpperBound(0) + 2); LowerUnBoundDimIndx[LowerUnBoundDimIndx.GetUpperBound(0)] = i; // MessageBox.Show("Dimension " + i + " is not lower-limited by any constraint." + "\n\r\n\r" + " Please check Problem: " + "\n\r\n\r" + OptimProblemTag + ".", // "MsgCode: 2019.04.02.A", MessageBoxButtons.OK, MessageBoxIcon.Information); // TimeToSolve_mSecs = Math.Round(Convert.ToDouble(DateTime.Now.ToString("ss.fffffff")) - TimeToSolve_mSecs, 5); } } } if (ProblemBoundnessCOndition) { // 3.3.2. Once identified an upper-constraint corresponding to a principal plane projection, // determine the closest lower-constraint. double[,] wProjFactr = new double[m, n]; int[][] RankedLowerConstrJndx = new int[n][]; double[][] LowerConstrDist = new double[n][]; string[][] OverlapCnstrSetW = new string[n][]; for (int i = 0; i <= n - 1; i++) { if (i != g) { Array.Resize(ref RankedLowerConstrJndx[i], 1); RankedLowerConstrJndx[i][0] = -1; Array.Resize(ref LowerConstrDist[i], 1); Array.Resize(ref OverlapCnstrSetW[i], 0); for (int j = 0; j <= m - 1; j++) { if (LowerBoundsDIMg[j, i]) { bool includeConstrWj = false; double wCritDistG = 0; int dimi = 0; // Detect Constraints perpendicular to axis dim i ... 2019.07.17 for (dimi = 0; dimi <= n - 1; dimi++){if (dimi != i) {if (A[j, dimi] != 0) { break;}}} if (dimi > n - 1) { // Constraint j is perpendicular to axis-dim i. Consider the value of dimension g where Constraint u cuts axis g wCritDistG = b[RankedUpperConstrJndx[i][0]] / A[RankedUpperConstrJndx[i][0], g]; includeConstrWj = true; } else { if (A[RankedUpperConstrJndx[i][0], i] == 0 || A[j, i] == 0 || A[j, g] == 0) { // Constraint is parallel to some axis. Either i or g if (A[j, g] == 0) { // wConstraint j is parallel to AxisDim g. Cuts on g at infinity // Constraint j is parallel to AxisDim g. Cut on i occurs at Gu = b[u] / A[u,i] ... 2019.07.17 // Conventional formula for wProjFactr[j, i] undefined for this case. Compute the angle formed by constraint j and principal plane ig // See and interprete drawing: in this case wShrtestDistG equals the value of coordinate Xcp[j, i]. wCritDistG = Xcp[j, i]; includeConstrWj = true; } else if (A[j, i] == 0) { // Projected wConstraint j is parallel to AxisDim i. Cuts on i at infinity wCritDistG = Xcp[j, g]; // Consider the value of dimension g where projected-Constraint j cuts axis g includeConstrWj = true; } else if (A[RankedUpperConstrJndx[i][0], i] == 0) { MessageBox.Show("Complete code here if necessary. 2019.04.06.A"); } else { MessageBox.Show("Why is it here?. 2019.07.18.A"); } } else { // Constraint is not parallel to any axis // As criterion use the distance d[j] for the COnstraint's j closest point Xcpj projected. See Figure 7. wCritDistG = Xcp[j, g] + Xcp[j, i] * A[j, g] / A[j, i]; includeConstrWj = true; } } if (includeConstrWj) { double AxisICutValue = b[j] / A[j, i]; double AxisGCutValue = b[j] / A[j, g]; int DegeneraIdx = 0; for (DegeneraIdx = 0; DegeneraIdx <= RankedLowerConstrJndx[i].GetUpperBound(0); DegeneraIdx++) { if (AxisICutValue == CutOverAxisI[i][DegeneraIdx] && AxisGCutValue == CutOverAxisG[i][DegeneraIdx]) { // Including this Constraint will cause Deneration of the Equation System. 2019.07.20 break; } } if (DegeneraIdx <= RankedLowerConstrJndx[i].GetUpperBound(0) || RankedLowerConstrJndx[i].GetUpperBound(0) < 0) { if (RankedLowerConstrJndx[i].GetUpperBound(0) < 0) { Array.Resize(ref OverlapCnstrSetW[i], OverlapCnstrSetW[i].GetUpperBound(0) + 2); OverlapCnstrSetW[i][OverlapCnstrSetW[i].GetUpperBound(0)] = j.ToString(); } else { OverlapCnstrSetW[i][DegeneraIdx] = OverlapCnstrSetW[i][DegeneraIdx] + "," + j.ToString(); } } else { Array.Resize(ref OverlapCnstrSetW[i], OverlapCnstrSetW[i].GetUpperBound(0) + 2); OverlapCnstrSetW[i][OverlapCnstrSetW[i].GetUpperBound(0)] = j.ToString(); } if (RankedLowerConstrJndx[i][0] >= 0) { Array.Resize(ref RankedLowerConstrJndx[i], RankedLowerConstrJndx[i].GetUpperBound(0) + 2); Array.Resize(ref LowerConstrDist[i], LowerConstrDist[i].GetUpperBound(0) + 2); } int AtJuIdx = 0; while (wCritDistG > LowerConstrDist[i][AtJuIdx] && AtJuIdx < RankedLowerConstrJndx[i].GetUpperBound(0)) { AtJuIdx++; } for (int SldIdx = RankedLowerConstrJndx[i].GetUpperBound(0); SldIdx > AtJuIdx; SldIdx--) { RankedLowerConstrJndx[i][SldIdx] = RankedLowerConstrJndx[i][SldIdx - 1]; LowerConstrDist[i][SldIdx] = LowerConstrDist[i][SldIdx - 1]; } RankedLowerConstrJndx[i][AtJuIdx] = j; LowerConstrDist[i][AtJuIdx] = wCritDistG; } } } } } // 3.4. Select active constraints. ------------------------------------------------------------------------------- int SelectIdx = -1; int[] CnstrSelctIdx = new int[0]; for (int i = 0; i <= n - 1; i++) { if (i != g) { // Verify overlaping projected-Constraints and avoid selecting Degenetative UpperConstraints int uprscn = 0; for (uprscn = 0; uprscn <= CnstrSelctIdx.GetUpperBound(0); uprscn++) { string[] OvrlpConstrSTRN = new string[1]; OvrlpConstrSTRN = OverlapCnstrSetU[i][0].Split(','); int[] OvrlpConstr = new int[OvrlpConstrSTRN.GetUpperBound(0) + 1]; for (int EvIdx = 0; EvIdx <= OvrlpConstrSTRN.GetUpperBound(0); EvIdx++) { OvrlpConstr[EvIdx] = Convert.ToInt32(OvrlpConstrSTRN[EvIdx]); } int EqvIdx = 0; for (EqvIdx = 0; EqvIdx <= OvrlpConstr.GetUpperBound(0); EqvIdx++) { int CnstrIdx = 0; for (CnstrIdx = 0; CnstrIdx <= CnstrSelctIdx.GetUpperBound(0); CnstrIdx++) { if(CnstrSelctIdx[CnstrIdx] == OvrlpConstr[EqvIdx]) { break; } } if(CnstrIdx <= CnstrSelctIdx.GetUpperBound(0)) { break; } } if (EqvIdx <= OvrlpConstr.GetUpperBound(0)) { break; } } if (uprscn > CnstrSelctIdx.GetUpperBound(0)) { Array.Resize(ref CnstrSelctIdx, CnstrSelctIdx.GetUpperBound(0) + 2); SelectIdx++; CnstrSelctIdx[SelectIdx] = RankedUpperConstrJndx[i][0]; } // Verify overlaping projected-Constraints and avoid selecting Degenetative LowerConstraints int wprscn = 0; for (wprscn = 0; wprscn <= OverlapCnstrSetW[i].GetUpperBound(0); wprscn++) { string[] OvrlpWConstrSTRN = new string[1]; OvrlpWConstrSTRN = OverlapCnstrSetW[i][wprscn].Split(','); int[] OvrlpWConstr = new int[OvrlpWConstrSTRN.GetUpperBound(0) + 1]; for (int EvIdx = 0; EvIdx <= OvrlpWConstrSTRN.GetUpperBound(0); EvIdx++) { OvrlpWConstr[EvIdx] = Convert.ToInt32(OvrlpWConstrSTRN[EvIdx]); } int EqvIdx = 0; for (EqvIdx = 0; EqvIdx <= OvrlpWConstr.GetUpperBound(0); EqvIdx++) { int CnstrIdx = 0; for (CnstrIdx = 0; CnstrIdx <= CnstrSelctIdx.GetUpperBound(0); CnstrIdx++) { if (CnstrSelctIdx[CnstrIdx] == OvrlpWConstr[EqvIdx]) { break; } } if (CnstrIdx <= CnstrSelctIdx.GetUpperBound(0)) { break; } } if (EqvIdx > OvrlpWConstr.GetUpperBound(0)) { Array.Resize(ref CnstrSelctIdx, CnstrSelctIdx.GetUpperBound(0) + 2); SelectIdx++; CnstrSelctIdx[SelectIdx] = OvrlpWConstr[wprscn]; break; } } } } ActiveKnstraints = CnstrSelctIdx[0].ToString(); for (int jSel = 1; jSel <= CnstrSelctIdx.GetUpperBound(0); jSel++) { ActiveKnstraints = ActiveKnstraints + "]0[" + CnstrSelctIdx[jSel].ToString(); } // 3.5 Find the extreme vertex’s coordinates by intersecting active constraints. --------------------------------- // Make BaseMatrixToInvert with Coeffs. Array.Resize(ref CnstrSelctIdx, n); double[][] BaseMatrixToInvert = TrisCSMatrixesCode.MatrixVers2.MatrixCreate(n, n); for (int jSel = 0; jSel <= CnstrSelctIdx.GetUpperBound(0); jSel++) { for (int iBase = 0; iBase <= n - 1; iBase++) BaseMatrixToInvert[jSel][iBase] = A[CnstrSelctIdx[jSel], iBase]; } // Matrix Inversion algorithm: MatrixInverse: https://jamesmccaffrey.wordpress.com/2015/03/06/inverting-a-matrix-using-c/ double[][] TheBaseMtrxInverse = TrisCSMatrixesCode.MatrixVers2.MatrixCreate(n, n); TheBaseMtrxInverse = TrisCSMatrixesCode.MatrixVers2.MatrixInverse(BaseMatrixToInvert); // Resource Vector. double[] TheResourceVector = TrisCSMatrixesCode.MatrixVers2.VectorCreate(n); for (int jSel = 0; jSel <= CnstrSelctIdx.GetUpperBound(0); jSel++) { TheResourceVector[jSel] = b[CnstrSelctIdx[jSel]]; } // Optimal Point. Constraints Intersection double[] OptimalPointCoords = TrisCSMatrixesCode.MatrixVers2.VectorCreate(n); OptimalPointCoords = TrisCSMatrixesCode.MatrixVers2.MatrixVectorProduct(TheBaseMtrxInverse, TheResourceVector); Array.Resize(ref x, n); for (int i = 0; i <= n - 1; i++) x[i] = OptimalPointCoords[i]; TimeToSolve_mSecs = Math.Round(Convert.ToDouble(DateTime.Now.ToString("ss.fffffff")) - TimeToSolve_mSecs, 5); // INI Round Values to Specified Precision int Precsn = 5; for (int i = 0; i <= C.GetUpperBound(0); i++) { C[i] = Math.Round(C[i], Precsn); x[i] = Math.Round(x[i], Precsn); } // INI Other Result Parameters // DominantDimension = g.ToString(); // Computed above. // UActive Constraints. Active Knstraints were Computed above. // Optimal Coordinates double OptimalValueFLOT = C[0] * x[0]; for (int i = 1; i <= n - 1; i++) { OptimalValueFLOT = OptimalValueFLOT + C[i] * x[i]; } OptimalValue = OptimalValueFLOT.ToString(); // bool[] InwardKnstraintIsRedundant = new bool[m]; // Inward Reduntant Constraints // bool[] OutwardKnstraintIsRedundant = new bool[m]; // Outward Reduntant Constraints } // Lower-unbounded dimensions if (LowerUnBoundDimIndx.GetUpperBound(0) >= 0){ LowerUnBoundDim = LowerUnBoundDimIndx[0].ToString(); for (int ubIdx = 1; ubIdx <= LowerUnBoundDimIndx.GetUpperBound(0); ubIdx++){ LowerUnBoundDim = LowerUnBoundDim + "]0[" + LowerUnBoundDimIndx[ubIdx].ToString(); // MessageBox.Show("Dimension " + i + " is not lower-limited by any constraint." + "\n\r\n\r" + " Please check Problem: " + "\n\r\n\r" + OptimProblemTag + ".", // "MsgCode: 2019.04.02.A", MessageBoxButtons.OK, MessageBoxIcon.Information); } TimeToSolve_mSecs = Math.Round(Convert.ToDouble(DateTime.Now.ToString("ss.fffffff")) - TimeToSolve_mSecs, 5); } } // Upper-unbounded dimensions if (UpperUnBoundDimIndx.GetUpperBound(0) >= 0) { UpperUnBoundDim = UpperUnBoundDimIndx[0].ToString(); for (int ubIdx = 1; ubIdx <= UpperUnBoundDimIndx.GetUpperBound(0); ubIdx++) { UpperUnBoundDim = UpperUnBoundDim + "]0[" + UpperUnBoundDimIndx[ubIdx].ToString(); // MessageBox.Show("Dimension " + i + " is not upper-limited by any constraint." + "\n\r\n\r" + " Please check Problem: " + "\n\r\n\r" + OptimProblemTag + ".", // "MsgCode: 2019.04.02.A", MessageBoxButtons.OK, MessageBoxIcon.Information); } TimeToSolve_mSecs = Math.Round(Convert.ToDouble(DateTime.Now.ToString("ss.fffffff")) - TimeToSolve_mSecs, 5); } } } }