1 /* 2 * ECOS - Embedded Conic Solver. 3 * Copyright (C) 2012-2015 A. Domahidi [domahidi@embotech.com], 4 * Automatic Control Lab, ETH Zurich & embotech GmbH, Zurich, Switzerland. 5 * 6 * This program is free software: you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation, either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with this program. If not, see <http://www.gnu.org/licenses/>. 18 */ 19 20 /* 21 * The exponential cone module is (c) Santiago Akle, Stanford University, 22 * [akle@stanford.edu] 23 */ 24 25 #include "expcone.h" 26 27 #ifdef EXPCONE 28 /* Evaluates the Hessian of the exponential dual cone barrier at the triplet 29 * w[0],w[1],w[2], and stores the upper triangular part of the matrix mu*H(w) 30 * at v[0],...,v[5] where the entries are arranged columnwise max_rows(pfloat * E,const spmat * mat)31 */ 32 void evalExpHessian(pfloat* w, pfloat* v, pfloat mu) 33 { 34 /** 35 * l = log(-y/x); 36 * r = -x*l-x+w; 37 * He = [[ 1/x^2 - 1/(r*x) + l^2/r^2, 1/(r*y) + (l*x)/(r^2*y), -l/r^2]; 38 * [ 1/(r*y) + (l*x)/(r^2*y), 1/y^2 - x/(r*y^2) + x^2/(r^2*y^2), -x/(r^2*y)]; 39 * [ -l/r^2, -x/(r^2*y), 1/r^2]]; 40 */ 41 42 pfloat x = w[0]; 43 pfloat y = w[1]; 44 pfloat z = w[2]; 45 pfloat l = log(-y/x); 46 pfloat r = -x*l-x+z; 47 v[0] = mu*((r*r-x*r+l*l*x*x)/(r*x*x*r)); 48 v[1] = mu*((z-x)/(r*r*y)); 49 v[2] = mu*((r*r-x*r+x*x)/(r*r*y*y)); 50 v[3] = mu*(-l/(r*r)); 51 v[4] = mu*(-x/(r*r*y)); 52 v[5] = mu*(1/(r*r)); 53 } sum_sq_rows(pfloat * E,const spmat * mat)54 55 /* Evaluates the gradient of the exponential cone g(z) at the triplet 56 * w[0],w[1],w[2], and stores the result at g[0],..,g[2] 57 */ 58 void evalExpGradient(pfloat* w, pfloat* g) 59 { 60 pfloat x = w[0]; 61 pfloat y = w[1]; 62 pfloat z = w[2]; 63 pfloat l = log(-y/x); 64 pfloat r = -x*l-x+z; 65 66 g[0] = (l*x-r)/(r*x); 67 g[1] = (x-r)/(y*r); 68 g[2] = -1/r; 69 } 70 71 /* Computes f_e(s_e) + f^\star_e(z_e) */ 72 pfloat evalBarrierValue(pfloat* siter, pfloat *ziter, idxint fc, idxint nexc) 73 { 74 pfloat l, u, v, w, x, y, z, o; 75 76 pfloat primal_barrier = 0.0; equilibrate_rows(const pfloat * E,spmat * mat)77 pfloat dual_barrier = 0.0; 78 79 idxint j; 80 81 /* Move to the first exponential cone slack */ 82 ziter = ziter+fc; 83 siter = siter+fc; 84 85 /* For the dual cone measure -u,v, -ul-u+w */ 86 /* For the primal cone measure z,v,omega-1 */ 87 for(j=0;j<nexc;j++) 88 { 89 /* Extract the entries */ equilibrate_cols(const pfloat * E,spmat * mat)90 u = ziter[0]; 91 v = ziter[1]; 92 w = ziter[2]; 93 94 x = siter[0]; 95 y = siter[1]; 96 z = siter[2]; 97 98 l = log(-v/u); 99 dual_barrier += -log(w-u-u*l)-log(-u)-log(v); 100 101 /* Primal Cone */ restore(const pfloat * D,const pfloat * E,spmat * mat)102 o = wrightOmega(1-x/z-log(z)+log(y)); 103 o = (o-1)*(o-1)/o; 104 primal_barrier += -log(o)-2*log(z)-log(y)-3; 105 106 ziter += 3; 107 siter += 3; 108 109 } 110 return primal_barrier+dual_barrier; 111 } 112 113 114 /* use_alternating_norm_equilibration(pwork * w)115 * Computes y[fc:end] += muH(x[fc:end])*x[fc:end], where 116 * fc is the index of the first exponential slack. 117 * This method assumes that the scalings have been updated by update scalings 118 * and that C->expc[cone_number].v contains mu*H(x). 119 * 120 */ 121 void scaleToAddExpcone(pfloat* y, pfloat* x, expcone* expc, idxint nexc, idxint fc) 122 { 123 idxint l; 124 /* Shift to the exponential slacks */ 125 x = x+fc; 126 y = y+fc; 127 128 for( l=0; l < nexc; l++ ){ 129 130 y[0]+= expc[l].v[0]*x[0]+expc[l].v[1]*x[1]+expc[l].v[3]*x[2]; 131 y[1]+= expc[l].v[1]*x[0]+expc[l].v[2]*x[1]+expc[l].v[4]*x[2]; 132 y[2]+= expc[l].v[3]*x[0]+expc[l].v[4]*x[1]+expc[l].v[5]*x[2]; 133 134 /* prepare index for next cone */ 135 x += 3; 136 y += 3; 137 } 138 } 139 140 /* 141 * Returns 1 if s is primal feasible 142 * with respect to the exponential cone, 143 * and 0 i.o.c 144 */ 145 idxint evalExpPrimalFeas(pfloat *s, idxint nexc) 146 { 147 pfloat x1,x2,x3,tmp1,psi; 148 idxint j = 0; 149 150 for(j =0 ; j < nexc; j++) 151 { 152 x1 = s[3*j]; 153 x2 = s[3*j+1]; 154 x3 = s[3*j+2]; 155 tmp1 = log(x2/x3); 156 psi = x3*tmp1 - x1; 157 if(psi<0||x2<0||x3<0) 158 { 159 return 0; 160 } 161 162 } 163 return 1; 164 } 165 166 /* 167 * Returns 1 if s is dual feasible 168 * with respect to the dual of the exponential cone, 169 * and 0 i.o.c 170 */ 171 idxint evalExpDualFeas(pfloat *z, idxint nexc) 172 { 173 pfloat x1,x2,x3,tmp1,psi; 174 idxint j = 0; 175 176 for(j =0 ; j < nexc; j++) 177 { 178 x1 = z[3*j]; 179 x2 = z[3*j+1]; 180 x3 = z[3*j+2]; 181 tmp1 = log(-x2/x1); 182 psi = -x1-x1*tmp1+x3; 183 if(0<x1||x2<0||psi<0) 184 { 185 return 0; 186 } 187 188 } 189 return 1; 190 } 191 192 #endif 193 194